.gitignore

@@ -8,4 +8,3 @@
 /.vscode
 /shared-folder
 /shared-folder/**
-/src/main/c/lib/**

.gitmodules

@@ -1,3 +0,0 @@
-[submodule "src/main/c/lib/stb"]
-	path = src/main/c/lib/stb
-	url = https://github.com/nothings/stb.git

Jenkinsfile

@@ -1,17 +0,0 @@
-pipeline {
-    agent any
-    tools {
-        maven '3.9.6'
-    }
-    stages {
-        stage('Build') {
-            steps {
-                sh 'mvn --version'
-                sh 'java -version'
-                sh 'apt-get update'
-                sh 'apt-get install gcc'
-                sh 'mvn clean generate-resources package'
-            }
-        }
-    }
-}

src/main/c/chunkmask.c

@@ -1,8 +1,7 @@
 #include <jni.h>
 #include <stdint.h>
-#include "../includes/libfluidsim.h"
+#include "includes/utilities.h"
-#include "../includes/utilities.h"
+#include "includes/chunkmask.h"
-#include "../includes/chunkmask.h"
 
 uint32_t matrix_transform(JNIEnv * env, jobjectArray jrx);
 
@@ -61,43 +60,43 @@ const uint32_t CHUNK_INDEX_ARR[] = {
 
 //control offsetting the advect sampler location if a valid neighbor chunk is hit
 const char CHUNK_NORMALIZE_U[] = {
-    1, 0, -1,
+    -1, 0, 1,
-    1, 0, -1,
+    -1, 0, 1,
-    1, 0, -1,
+    -1, 0, 1,
 
-    1, 0, -1,
+    -1, 0, 1,
-    1, 0, -1,
+    -1, 0, 1,
-    1, 0, -1,
+    -1, 0, 1,
 
-    1, 0, -1,
+    -1, 0, 1,
-    1, 0, -1,
+    -1, 0, 1,
-    1, 0, -1,
+    -1, 0, 1,
 };
 
 const char CHUNK_NORMALIZE_V[] = {
-     1,  1,  1,
-     0,  0,  0,
     -1, -1, -1,
+     0,  0,  0,
+     1,  1,  1,
 
-    1,  1,  1,
-     0,  0,  0,
     -1, -1, -1,
+     0,  0,  0,
+     1,  1,  1,
 
-    1,  1,  1,
-     0,  0,  0,
     -1, -1, -1,
+     0,  0,  0,
+     1,  1,  1,
 };
 
 const char CHUNK_NORMALIZE_W[] = {
-    1, 1, 1,
+    -1, -1, -1,
-    1, 1, 1,
+    -1, -1, -1,
-    1, 1, 1,
+    -1, -1, -1,
 
     0, 0, 0,
     0, 0, 0,
     0, 0, 0,
 
-    -1, -1, -1,
+    1, 1, 1,
-    -1, -1, -1,
+    1, 1, 1,
-    -1, -1, -1,
+    1, 1, 1,
 };

src/main/c/compile.sh

@@ -42,17 +42,22 @@ rm -f ./*.dll
 
 #compile object files
 COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
-INPUT_FILES="./src/densitystep.c"
+INPUT_FILES="./fluidsim.c"
+OUTPUT_FILE="./fluidsim.o"
+gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
+
+COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
+INPUT_FILES="./densitystep.c"
 OUTPUT_FILE="./densitystep.o"
 gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
 
 COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
-INPUT_FILES="./src/velocitystep.c"
+INPUT_FILES="./velocitystep.c"
 OUTPUT_FILE="./velocitystep.o"
 gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
 
 COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
-INPUT_FILES="./src/chunkmask.c"
+INPUT_FILES="./chunkmask.c"
 OUTPUT_FILE="./chunkmask.o"
 gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
 

src/main/c/densitystep.c

@@ -2,13 +2,11 @@
 #include <stdio.h>
 #include <immintrin.h>
 #include <stdint.h>
-#include <pthread.h>
 
-#include "../includes/libfluidsim.h"
+#include "includes/utilities.h"
-#include "../includes/utilities.h"
+#include "includes/chunkmask.h"
-#include "../includes/chunkmask.h"
 
-void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, jobjectArray d0, float * u, float * v, float * w, float 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);
 
 /*
  * Class:     electrosphere_FluidSim
@@ -21,12 +19,12 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_addDensity
   jint N,
   jint chunk_mask,
   jobjectArray jrx,
-  jobjectArray x0,
+  jobject x0,
   jfloat dt){
     int i;
     int size=N*N*N;
     float * x = GET_ARR(env,jrx,CENTER_LOC);
-    float * s = GET_ARR(env,x0,CENTER_LOC);
+    float * s = (*env)->GetDirectBufferAddress(env,x0);
     for(i=0; i<size; i++){
         x[i] += dt*s[i];
     }
@@ -43,7 +41,7 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
   jint N,
   jint chunk_mask,
   jobjectArray jrx,
-  jobjectArray jrx0,
+  jobject jrx0,
   jobjectArray jru,
   jobjectArray jrv,
   jobjectArray jrw,
@@ -54,7 +52,7 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
     float c=1+6*a;
     int i, j, k, l, m;
     float * x = GET_ARR(env,jrx,CENTER_LOC);
-    float * x0 = GET_ARR(env,jrx0,CENTER_LOC);
+    float * x0 = (*env)->GetDirectBufferAddress(env,jrx0);
     
     __m256 aScalar = _mm256_set1_ps(a);
     __m256 cScalar = _mm256_set1_ps(c);
@@ -97,17 +95,18 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectDensity
   jint N,
   jint chunk_mask,
   jobjectArray jrx,
-  jobjectArray jrx0,
+  jobject jrx0,
   jobjectArray jru,
   jobjectArray jrv,
   jobjectArray jrw,
   jfloat DIFFUSION_CONST,
   jfloat VISCOSITY_CONST,
   jfloat dt){
-    advectDensity(env,chunk_mask,N,3,jrx,jrx0,GET_ARR(env,jru,CENTER_LOC),GET_ARR(env,jrv,CENTER_LOC),GET_ARR(env,jrw,CENTER_LOC),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, jobjectArray jrd0, float * u, float * v, float * w, float 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;
@@ -116,12 +115,12 @@ void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray
 
     float * d = GET_ARR(env,jrd,CENTER_LOC);
 
-    float * d0 = GET_ARR(env,jrd0,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++){
-                d0 = GET_ARR(env,jrd0,CENTER_LOC);
+                sampleArr = d0;
                 //calculate location to pull from
                 x = i-dtx*u[IX(i,j,k)];
                 y = j-dty*v[IX(i,j,k)];
@@ -129,12 +128,12 @@ void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray
 
                 m = n = o = 1;
 
-                // if(x < 1){ m -= 1; }
+                if(x < 1){ m -= 1; }
-                // if(x >= N-1){ m += 1; }
+                if(x >= N-1){ m += 1; }
-                // if(y < 1){ n -= 1; }
+                if(y < 1){ n -= 1; }
-                // if(y >= N-1){ n += 1; }
+                if(y >= N-1){ n += 1; }
-                // if(z < 1){ o -= 1; }
+                if(z < 1){ o -= 1; }
-                // if(z >= N-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){
@@ -145,56 +144,31 @@ void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray
                 // }
                 // if(CK(m,n,o) != CENTER_LOC && ARR_EXISTS(chunk_mask,m,n,o)){
                 //     // printf("Hit other chunk\n");
-                //     d0 = GET_ARR(env,jrd0,CK(m,n,o));
+                //     sampleArr = GET_ARR(env,jrd,CK(m,n,o));
-                //     x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * (N-1);
+                //     x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * N;
-                //     y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * (N-1);
+                //     y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * N;
-                //     z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * (N-1);
+                //     z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * N;
                 // }
 
-                if(x < 0.001f){
-                    //cases to consider:
-                    //m = 0, x = -10
-                    //m = 2, x = 0.01
-                    x=0.001f;
-                    i0=(int)0;
-                    i1=1;
-                    s0 = 0.999f;
-                    s1 = 0.001f;
-                } else if(x >= N - 1){
-                    //cases to consider:
-                    //m = 0, x = 17.01
-                    //m = 2, x = 20
-                    x = N-1;
-                    i0=(int)N-2;
-                    i1=N-1;
-                    s0 = 0.001f;
-                    s1 = 0.999f;
-                } else {
-                    i0=(int)x;
-                    i1=i0+1;
-                    s1 = x-i0;
-                    s0 = 1-s1;
-                }
-
                 //clamp location within chunk
-                // if (x<0.5f) x=0.5f;
+                if (x<0.5f) x=0.5f;
-                // if (x>N+0.5f) x=N+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;
+                i0=(int)x;
-                // i1=i0+1;
+                i1=i0+1;
                 j0=(int)y;
                 j1=j0+1;
                 k0=(int)z;
                 k1=k0+1;
 
                 //calculate percentage of each index
-                // s1 = x-i0;
+                s1 = x-i0;
-                // s0 = 1-s1;
+                s0 = 1-s1;
                 t1 = y-j0;
                 t0 = 1-t1;
                 u1 = z-k0;
@@ -238,16 +212,16 @@ void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray
                 // }
                 d[IX(i,j,k)] = 
                 s0*(
-                    t0*u0*d0[IX(i0,j0,k0)]+
+                    t0*u0*sampleArr[IX(i0,j0,k0)]+
-                    t1*u0*d0[IX(i0,j1,k0)]+
+                    t1*u0*sampleArr[IX(i0,j1,k0)]+
-                    t0*u1*d0[IX(i0,j0,k1)]+
+                    t0*u1*sampleArr[IX(i0,j0,k1)]+
-                    t1*u1*d0[IX(i0,j1,k1)]
+                    t1*u1*sampleArr[IX(i0,j1,k1)]
                 )+
                 s1*(
-                    t0*u0*d0[IX(i1,j0,k0)]+
+                    t0*u0*sampleArr[IX(i1,j0,k0)]+
-                    t1*u0*d0[IX(i1,j1,k0)]+
+                    t1*u0*sampleArr[IX(i1,j1,k0)]+
-                    t0*u1*d0[IX(i1,j0,k1)]+
+                    t0*u1*sampleArr[IX(i1,j0,k1)]+
-                    t1*u1*d0[IX(i1,j1,k1)]
+                    t1*u1*sampleArr[IX(i1,j1,k1)]
                 );
             }
         }

src/main/c/fluidsim.c

@@ -0,0 +1,396 @@
+#include <jni.h>
+#include <stdio.h>
+#include <immintrin.h>
+#include <stdint.h>
+
+#include "includes/utilities.h"
+#include "includes/chunkmask.h"
+
+
+
+#define LINEARSOLVERTIMES 10
+
+
+
+void diffuse(JNIEnv * env, uint32_t chunk_mask, int N, int b, float * x, float * x0, float diff, 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);
+void project(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float * p, float * div);
+void set_bnd(JNIEnv * env, uint32_t chunk_mask, int N, int b, float * x);
+void dens_step(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jrx, float * x0, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float diff, float dt);
+void vel_step(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float * u0, float * v0, float * w0, float visc, float dt);
+void lin_solve(JNIEnv * env, uint32_t chunk_mask, int N, int b, float* x, float* x0, float a, float c);
+
+/**
+ * The core simulation function
+*/
+JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
+    JNIEnv * env,
+    jobject this,
+    jint DIM_X,
+    jint chunk_mask,
+    jobject jx,
+    jobject jx0,
+    jobject ju,
+    jobject jv,
+    jobject jw,
+    jobject ju0,
+    jobject jv0,
+    jobject jw0,
+    jfloat DIFFUSION_RATE,
+    jfloat VISCOSITY_RATE,
+    jfloat timestep){
+    jboolean isCopy;
+    // float * x = (*env)->GetDirectBufferAddress(env,jx);
+    float * x0 = (*env)->GetDirectBufferAddress(env,jx0);
+    // float * u = (*env)->GetDirectBufferAddress(env,ju);
+    // float * v = (*env)->GetDirectBufferAddress(env,jv);
+    // float * w = (*env)->GetDirectBufferAddress(env,jw);
+    float * u0 = (*env)->GetDirectBufferAddress(env,ju0);
+    float * v0 = (*env)->GetDirectBufferAddress(env,jv0);
+    float * w0 = (*env)->GetDirectBufferAddress(env,jw0);
+    int N = DIM_X;
+    int i,j,k;
+    vel_step(env, chunk_mask, DIM_X, ju, jv, jw, u0, v0, w0, VISCOSITY_RATE, timestep);
+    dens_step(env, chunk_mask, DIM_X, jx, x0, ju, jv, jw, DIFFUSION_RATE, timestep);
+}
+
+/**
+ * Adds values from a source array to a current frame array (eg more density to the main density array)
+*/
+void add_source(int N, float * x, float * s, float dt){
+	int i;
+    int size=N*N*N;
+	for(i=0; i<size; i++){
+        x[i] += dt*s[i];
+    }
+}
+
+/**
+ * Diffuses a given array by a diffusion constant
+*/
+void diffuse(JNIEnv * env, uint32_t chunk_mask, int N, int b, float * x, float * x0, float diff, float dt){
+    float a=dt*diff*N*N*N;
+    lin_solve(env, chunk_mask, N, b, x, x0, a, 1+6*a);
+}
+
+/**
+ * Advects a given array based on the force vectors in the simulation
+*/
+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;
+    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);
+
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            for(i=1; i<N-1; i++){
+                x = i-dtx*u[IX(i,j,k)]; y = j-dty*v[IX(i,j,k)]; z = k-dtz*w[IX(i,j,k)];
+                if (x<0.5f) x=0.5f; if (x>N+0.5f) x=N+0.5f; i0=(int)x; i1=i0+1;
+                if (y<0.5f) y=0.5f; if (y>N+0.5f) y=N+0.5f; j0=(int)y; j1=j0+1;
+                if (z<0.5f) z=0.5f; if (z>N+0.5f) z=N+0.5f; k0=(int)z; k1=k0+1;
+
+                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*d0[IX(i0,j0,k0)]+t1*u0*d0[IX(i0,j1,k0)]+t0*u1*d0[IX(i0,j0,k1)]+t1*u1*d0[IX(i0,j1,k1)])+
+                    s1*(t0*u0*d0[IX(i1,j0,k0)]+t1*u0*d0[IX(i1,j1,k0)]+t0*u1*d0[IX(i1,j0,k1)]+t1*u1*d0[IX(i1,j1,k1)]);
+            }
+        }
+    }
+    set_bnd(env, chunk_mask, N, b, d);
+}
+
+/**
+ * The main density step function
+*/
+void dens_step(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jrx, float * x0, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float diff, float dt){
+    float * x = GET_ARR(env,jrx,CENTER_LOC);
+    float * u = GET_ARR(env,jru,CENTER_LOC);
+    float * v = GET_ARR(env,jrv,CENTER_LOC);
+    float * w = GET_ARR(env,jrw,CENTER_LOC);
+    add_source(N, x, x0, dt);
+    SWAP(x0, x);
+    diffuse(env, chunk_mask, N, 0, x, x0, diff, dt);
+    SWAP(x0, x);
+    advect(env, chunk_mask, N, 0, jrx, x0, u, v, w, dt);
+}
+
+/**
+ * The main velocity step function
+*/
+void vel_step(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float * u0, float * v0, float * w0, float visc, float dt){
+    float * u = GET_ARR(env,jru,CENTER_LOC);
+    float * v = GET_ARR(env,jrv,CENTER_LOC);
+    float * w = GET_ARR(env,jrw,CENTER_LOC);
+    add_source(N, u, u0, dt);
+    add_source(N, v, v0, dt);
+    add_source(N, w, w0, dt);
+    SWAP(u0, u);
+    diffuse(env, chunk_mask, N, 1, u, u0, visc, dt);
+    SWAP(v0, v);
+    diffuse(env, chunk_mask, N, 2, v, v0, visc, dt);
+    SWAP(w0, w);
+    diffuse(env, chunk_mask, N, 3, w, w0, visc, dt);
+    project(env, chunk_mask, N, jru, jrv, jrw, u0, v0);
+    SWAP(u0, u);
+    SWAP(v0, v);
+    SWAP(w0, w);
+    advect(env, chunk_mask, N, 1, jru, u0, u0, v0, w0, dt);
+    advect(env, chunk_mask, N, 2, jrv, v0, u0, v0, w0, dt);
+    advect(env, chunk_mask, N, 3, jrw, w0, u0, v0, w0, dt);
+    project(env, chunk_mask, N, jru, jrv, jrw, u0, v0);
+}
+
+//used for temporary vector storage when appropriate
+float container[16];
+
+/**
+ * Projects a given array based on force vectors
+*/
+void project(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float * p, float * div){
+    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);
+
+    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, 0, div);
+    set_bnd(env, chunk_mask, N, 0, p);
+
+    lin_solve(env, chunk_mask, N, 0, p, div, 1, 6);
+
+
+    constScalar = _mm256_set1_ps(0.5f*N);
+    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(env, chunk_mask, N, 1, u);
+    set_bnd(env, chunk_mask, N, 2, v);
+    set_bnd(env, chunk_mask, N, 3, w);
+}
+
+/**
+ * Solves a linear system of equations in a vectorized manner
+*/
+void lin_solve(JNIEnv * env, uint32_t chunk_mask, int N, int b, float* x, float* x0, float a, float c){
+    int i, j, k, l, m;
+    __m256 aScalar = _mm256_set1_ps(a);
+    __m256 cScalar = _mm256_set1_ps(c);
+    // iterate the solver
+    for ( l=0 ; l<LINEARSOLVERTIMES ; l++ ) {
+        // 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, N, b, x);
+    }
+}
+
+/**
+ * Sets the bounds of the simulation
+*/
+void set_bnd(JNIEnv * env, uint32_t chunk_mask, int N, int b, float * target){
+    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/c/includes/electrosphere_FluidSim.h

@@ -10,13 +10,13 @@ extern "C" {
 #undef electrosphere_FluidSim_DIM
 #define electrosphere_FluidSim_DIM 18L
 #undef electrosphere_FluidSim_DIFFUSION_CONSTANT
-#define electrosphere_FluidSim_DIFFUSION_CONSTANT 0.0f
+#define electrosphere_FluidSim_DIFFUSION_CONSTANT 1.0E-5f
 #undef electrosphere_FluidSim_VISCOSITY_CONSTANT
-#define electrosphere_FluidSim_VISCOSITY_CONSTANT 0.0f
+#define electrosphere_FluidSim_VISCOSITY_CONSTANT 1.0E-5f
 #undef electrosphere_FluidSim_LINEARSOLVERTIMES
-#define electrosphere_FluidSim_LINEARSOLVERTIMES 20L
+#define electrosphere_FluidSim_LINEARSOLVERTIMES 10L
 #undef electrosphere_FluidSim_GRAVITY
-#define electrosphere_FluidSim_GRAVITY -100.0f
+#define electrosphere_FluidSim_GRAVITY -1000.0f
 /*
  * Class:     electrosphere_FluidSim
  * Method:    simulate
@@ -84,26 +84,26 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectVectors
 /*
  * Class:     electrosphere_FluidSim
  * Method:    addDensity
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;F)V
+ * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;F)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_addDensity
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
  * Method:    solveDiffuseDensity
- * Signature: (II[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;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
  * Method:    advectDensity
- * Signature: (II[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;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectDensity
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
@@ -113,14 +113,6 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectDensity
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setBoundsToNeighbors
   (JNIEnv *, jobject, jint, jint, jint, jobjectArray);
 
-/*
- * Class:     electrosphere_FluidSim
- * Method:    copyNeighbors
- * Signature: (IIII[Ljava/nio/ByteBuffer;)V
- */
-JNIEXPORT void JNICALL Java_electrosphere_FluidSim_copyNeighbors
-  (JNIEnv *, jobject, jint, jint, jint, jint, jobjectArray);
-
 #ifdef __cplusplus
 }
 #endif

src/main/c/includes/libfluidsim.h

@@ -1,10 +0,0 @@
-//include guard
-#ifndef LIB_FLUID_SIM
-#define LIB_FLUID_SIM
-
-//include stb ds
-#define STB_IMAGE_IMPLEMENTATION
-#include "../lib/stb/stb_ds.h"
-
-//close include guard
-#endif

src/main/c/lib/stb

@@ -1 +0,0 @@
-Subproject commit ae721c50eaf761660b4f90cc590453cdb0c2acd0

src/main/c/velocitystep.c

@@ -3,9 +3,8 @@
 #include <immintrin.h>
 #include <stdint.h>
 
-#include "../includes/libfluidsim.h"
+#include "includes/utilities.h"
-#include "../includes/utilities.h"
+#include "includes/chunkmask.h"
-#include "../includes/chunkmask.h"
 
 
 #define BOUND_NO_DIR 0
@@ -17,7 +16,7 @@
 #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, jobjectArray jrd0, float * u, float * v, float * w, 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);
 
 
 /*
@@ -67,7 +66,7 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveVectorDiffuse
   jfloat DIFFUSION_CONST,
   jfloat VISCOSITY_CONST,
   jfloat dt){
-    float a=dt*VISCOSITY_CONST*N*N*N;
+    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);
@@ -178,9 +177,7 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupProjection
   jfloat dt){
     int i, j, k;
 
-    __m256 xVector = _mm256_set1_ps(N);
+    __m256 nVector = _mm256_set1_ps(N);
-    __m256 yVector = _mm256_set1_ps(N);
-    __m256 zVector = _mm256_set1_ps(N);
     __m256 constScalar = _mm256_set1_ps(-1.0/3.0);
     __m256 zeroVec = _mm256_set1_ps(0);
     __m256 vector, vector2, vector3;
@@ -192,9 +189,6 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupProjection
     float * p = GET_ARR(env,jru0,CENTER_LOC);
     float * div = GET_ARR(env,jrv0,CENTER_LOC);
 
-    float scalar = 1.0/3.0;
-    float h = 1.0/N;
-
     for(k=1; k<N-1; k++){
         for(j=1; j<N-1; j++){
             i = 1;
@@ -204,15 +198,15 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupProjection
             //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,xVector);
+            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,yVector);
+            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,zVector);
+            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);
@@ -226,30 +220,21 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupProjection
             //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,xVector);
+            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,yVector);
+            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,zVector);
+            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);
-            
-
-            // for(i = 1; i < N - 1; i++){
-            //     div[IX(i,j,k)] = 
-            //     -scalar*h*(u[IX(i+1,j,k)]-u[IX(i-1,j,k)]+
-            //     v[IX(i,j+1,k)]-v[IX(i,j-1,k)]+
-            //     w[IX(i,j,k+1)]-w[IX(i,j,k-1)]);
-            //     p[IX(i,j,k)] = 0;
-            // }
         }
     }
 }
@@ -277,34 +262,31 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveProjection
     __m256 aScalar = _mm256_set1_ps(a);
     __m256 cScalar = _mm256_set1_ps(c);
 
-    float * p = GET_ARR(env,jru0,CENTER_LOC);
+    float * x = GET_ARR(env,jru0,CENTER_LOC);
-    float * div = GET_ARR(env,jrv0,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(&p[IX(i-1,j,k)]);
+                __m256 vector = _mm256_loadu_ps(&x[IX(i-1,j,k)]);
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[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(&p[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(&p[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(&p[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k-1)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[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_mul_ps(vector,aScalar);
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&div[IX(i,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x0[IX(i,j,k)]));
                 vector = _mm256_div_ps(vector,cScalar);
-                _mm256_storeu_ps(&p[IX(i,j,k)],vector);
+                _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++){
-                    p[IX(i,j,k)] = (div[IX(i,j,k)] + a*(p[IX(i-1,j,k)]+p[IX(i+1,j,k)]+p[IX(i,j-1,k)]+p[IX(i,j+1,k)]+p[IX(i,j,k-1)]+p[IX(i,j,k+1)]))/c;
+                    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;
                 }
             }
-            // for(i=1; i < N-1; i++){
-            //     p[IX(i,j,k)] = (div[IX(i,j,k)] + a*(p[IX(i-1,j,k)]+p[IX(i+1,j,k)]+p[IX(i,j-1,k)]+p[IX(i,j+1,k)]+p[IX(i,j,k-1)]+p[IX(i,j,k+1)]))/c;
-            // }
         }
     }
 }
@@ -327,10 +309,9 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection
   jfloat VISCOSITY_CONST,
   jfloat dt){
     int i, j, k;
-    // __m256 constScalar = _mm256_set1_ps(0.5f*N);
+    __m256 nVector = _mm256_set1_ps(N);
-    __m256 xScalar = _mm256_set1_ps(0.5*N);
+    __m256 constScalar = _mm256_set1_ps(0.5f*N);
-    __m256 yScalar = _mm256_set1_ps(0.5*N);
+    __m256 zeroVec = _mm256_set1_ps(0);
-    __m256 zScalar = _mm256_set1_ps(0.5*N);
     __m256 vector, vector2, vector3;
 
     float * u = GET_ARR(env,jru,CENTER_LOC);
@@ -340,8 +321,6 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection
     float * p = GET_ARR(env,jru0,CENTER_LOC);
     float * div = GET_ARR(env,jrv0,CENTER_LOC);
 
-    float h = 1.0 / N;
-
     for ( k=1 ; k<N-1 ; k++ ) {
         for ( j=1 ; j<N-1 ; j++ ) {
             //
@@ -351,14 +330,14 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection
             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,xScalar);
+            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,xScalar);
+            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);
             //
@@ -368,14 +347,14 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection
             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,yScalar);
+            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,yScalar);
+            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);
             //
@@ -385,21 +364,16 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection
             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,zScalar);
+            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,zScalar);
+            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);
-            // for(i = 1; i < N-1; i++){
-            //     u[IX(i,j,k)] = u[IX(i,j,k)] - 0.5 * (p[IX(i+1,j,k)] - p[IX(i-1,j,k)]) / h;
-            //     v[IX(i,j,k)] = v[IX(i,j,k)] - 0.5 * (p[IX(i,j+1,k)] - p[IX(i,j-1,k)]) / h;
-            //     w[IX(i,j,k)] = w[IX(i,j,k)] - 0.5 * (p[IX(i,j,k+1)] - p[IX(i,j,k-1)]) / h;
-            // }
         }
     }
 }
@@ -421,13 +395,13 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectVectors
   jfloat DIFFUSION_CONST,
   jfloat VISCOSITY_CONST,
   jfloat dt){
-    advect(env,chunk_mask,N,1,jru,jru0,GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),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,jrv0,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,jrw0,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, jobjectArray jrd0, float * u, float * v, float * w, 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){
     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;
@@ -436,12 +410,12 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, j
 
     float * d = GET_ARR(env,jrd,CENTER_LOC);
 
-    float * d0 = GET_ARR(env,jrd0,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++){
-                d0 = GET_ARR(env,jrd0,CENTER_LOC);
+                sampleArr = d0;
                 //calculate location to pull from
                 x = i-dtx*u[IX(i,j,k)];
                 y = j-dty*v[IX(i,j,k)];
@@ -449,175 +423,55 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, j
 
                 m = n = o = 1;
 
-                if(x < 0){ m += 1; }
+                if(x < 1){ m -= 1; }
-                else if(x >= N){ m -= 1; }
+                if(x >= N-1){ m += 1; }
-                if(y < 0){ n += 1; }
+                if(y < 1){ n -= 1; }
-                else if(y >= N){ n -= 1; }
+                if(y >= N-1){ n += 1; }
-                if(z < 0){ o += 1; }
+                if(z < 1){ o -= 1; }
-                else if(z >= N){ 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 && ARR_EXISTS(chunk_mask,m,n,o)){
+                // if(CK(m,n,o) != CENTER_LOC){
-
+                //     printf("Looking in border chunk\n");
-                    // if(i == 1 && j == 1 && k == 1){
+                // }
-                    //     printf("\narr indices: %d %d %d\n\n",m,n,o);
+                // 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;
                 // }
-
-                    //cases:
-                    //if x = 17.01, m = 2
-                    // 17 in current array is 1 in neighbor
-                    // 18 in current array is 2 in neighbor
-                    // 19 in current array is 3 in neighbor
-                    //want to sample neighbor array at 1 & 2
-                    //x becomes 1.01, sampling new array (keep in mind that 0 in the new array should contain the current array values)
-                    //modification: subtract 16
-
-                    //cases:
-                    //if x = 16.99, m = 2
-                    // 16 in current array is 0 in neighbor
-                    // 17 in current array is 1 in neighbor
-                    // 18 in current array is 2 in neighbor
-                    // 19 in current array is 3 in neighbor
-                    //want to sample current array still
-                    //x becomes 1.01, sampling new array (keep in mind that 0 in the new array should contain the current array values)
-                    //modification: no modification
-
-                    //if x = 0.01, m = 0
-                    // 0 in current array is 16 in neighbor
-                    //-1 in current array is 15 in neighbor
-                    //-2 in current array is 14 in neighbor
-                    //want to sample current array still
-                    //x becomes 15.01, sampling new array (keep in mind that 17 in the new array should contain the current array)
-                    //modification: no modification
-
-                    //if x = -0.01, m = 0
-                    // 0 in current array is 16 in neighbor
-                    //-1 in current array is 15 in neighbor
-                    //-2 in current array is 14 in neighbor
-                    //want to sample -1 & 0, so i0 becomes 15
-                    //x becomes 15.99, sampling new array (keep in mind that 17 in the new array should contain the current array)
-                    //modification: add 16
-
-                    //if x = -2, m = 0
-                    // 0 in current array is 16 in neighbor
-                    //-1 in current array is 15 in neighbor
-                    //-2 in current array is 14 in neighbor
-                    //x becomes 14, sampling new array (keep in mind that 17 in the new array should contain the current array)
-                    //modification: add 16
-
-
-                    // printf("Hit other chunk\n");
-                    d0 = GET_ARR(env,jrd0,CK(m,n,o));
-                    x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * (N-2);
-                    // printf("%d => %f\n",m,x);
-                    y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * (N-2);
-                    z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * (N-2);
-                }
 
                 //clamp location within chunk
-                //get indices, and calculate percentage to pull from each index
+                if (x<0.5f) x=0.5f;
-                if(x < 0.001f){
+                if (x>N+0.5f) x=N+0.5f;
-                    //cases to consider:
+                if (y<0.5f) y=0.5f;
-                    //m = 0, x = -10
+                if (y>N+0.5f) y=N+0.5f;
-                    //m = 2, x = 0.01
+                if (z<0.5f) z=0.5f;
-                    x=0.001f;
+                if (z>N+0.5f) z=N+0.5f;
-                    i0=(int)0;
-                    i1=1;
-                    s0 = 0.999f;
-                    s1 = 0.001f;
-                } else if(x > N - 1){
-                    //cases to consider:
-                    //m = 0, x = 17.01
-                    //m = 2, x = 20
-                    x = N-1;
-                    i0=(int)N-2;
-                    i1=N-1;
-                    s0 = 0.001f;
-                    s1 = 0.999f;
-                } else {
-                    i0=(int)x;
-                    i1=i0+1;
-                    s1 = x-i0;
-                    s0 = 1-s1;
-                }
-
-                if(y < 0.001f){
-                    //cases to consider:
-                    //m = 0, x = -10
-                    //m = 2, x = 0.01
-                    y=0.001f;
-                    j0=(int)0;
-                    j1=1;
-                    t0 = 0.999f;
-                    t1 = 0.001f;
-                } else if(y > N - 1){
-                    //cases to consider:
-                    //m = 0, x = 17.01
-                    //m = 2, x = 20
-                    y = N-1;
-                    j0=(int)N-2;
-                    j1=N-1;
-                    t0 = 0.001f;
-                    t1 = 0.999f;
-                } else {
-                    j0=(int)y;
-                    j1=j0+1;
-                    t1 = y-j0;
-                    t0 = 1-t1;
-                }
-
-
-                if(z < 0.001f){
-                    //cases to consider:
-                    //m = 0, x = -10
-                    //m = 2, x = 0.01
-                    z=0.001f;
-                    k0=(int)0;
-                    k1=1;
-                    u0 = 0.999f;
-                    u1 = 0.001f;
-                } else if(z > N - 1){
-                    //cases to consider:
-                    //m = 0, x = 17.01
-                    //m = 2, x = 20
-                    z = N-1;
-                    k0=(int)N-2;
-                    k1=N-1;
-                    u0 = 0.001f;
-                    u1 = 0.999f;
-                } else {
-                    k0=(int)z;
-                    k1=k0+1;
-                    u1 = z-k0;
-                    u0 = 1-u1;
-                }
-
-                // if (x<0.001f) x=0.001f;
-                // if (x>N+0.5f) x=N+0.5f;
-                // if (y<0.001f) y=0.001f;
-                // if (y>N+0.5f) y=N+0.5f;
-                // if (z<0.001f) z=0.001f;
-                // if (z>N+0.5f) z=N+0.5f;
 
                 //get actual indices
-                // i0=(int)x;
+                i0=(int)x;
-                // i1=i0+1;
+                i1=i0+1;
-                // j0=(int)y;
+                j0=(int)y;
-                // j1=j0+1;
+                j1=j0+1;
-                // k0=(int)z;
+                k0=(int)z;
-                // k1=k0+1;
+                k1=k0+1;
 
                 //calculate percentage of each index
-                // s1 = x-i0;
+                s1 = x-i0;
-                // s0 = 1-s1;
+                s0 = 1-s1;
-                // t1 = y-j0;
+                t1 = y-j0;
-                // t0 = 1-t1;
+                t0 = 1-t1;
-                // u1 = z-k0;
+                u1 = z-k0;
-                // u0 = 1-u1;
+                u0 = 1-u1;
 
-                // if(i0 >= N){
+                if(i0 >= N){
-                //     i0 = N - 1;
+                    i0 = N - 1;
-                // }
+                }
                 // if(i0 < 0){
                 //     i0 = 0;
                 // }
@@ -633,9 +487,9 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, j
                 // if(k0 < 0){
                 //     k0 = 0;
                 // }
-                // if(i1 >= N){
+                if(i1 >= N){
-                //     i1 = N - 1;
+                    i1 = N - 1;
-                // }
+                }
                 // if(i1 < 0){
                 //     i1 = 0;
                 // }
@@ -653,16 +507,16 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, j
                 // }
                 d[IX(i,j,k)] = 
                 s0*(
-                    t0*u0*d0[IX(i0,j0,k0)]+
+                    t0*u0*sampleArr[IX(i0,j0,k0)]+
-                    t1*u0*d0[IX(i0,j1,k0)]+
+                    t1*u0*sampleArr[IX(i0,j1,k0)]+
-                    t0*u1*d0[IX(i0,j0,k1)]+
+                    t0*u1*sampleArr[IX(i0,j0,k1)]+
-                    t1*u1*d0[IX(i0,j1,k1)]
+                    t1*u1*sampleArr[IX(i0,j1,k1)]
                 )+
                 s1*(
-                    t0*u0*d0[IX(i1,j0,k0)]+
+                    t0*u0*sampleArr[IX(i1,j0,k0)]+
-                    t1*u0*d0[IX(i1,j1,k0)]+
+                    t1*u0*sampleArr[IX(i1,j1,k0)]+
-                    t0*u1*d0[IX(i1,j0,k1)]+
+                    t0*u1*sampleArr[IX(i1,j0,k1)]+
-                    t1*u1*d0[IX(i1,j1,k1)]
+                    t1*u1*sampleArr[IX(i1,j1,k1)]
                 );
             }
         }
@@ -679,16 +533,35 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setBoundsToNeighbors
     int DIM = N;
     float * target = GET_ARR(env,neighborArray,CENTER_LOC);
     float * source;
+    // if(ARR_EXISTS(chunk_mask,0,1,1)){
+    //     source = GET_ARR(env,neighborArray,CK(0,1,1));
+    //     for(int x=1; x < DIM-1; x++){
+    //         for(int y = 1; y < DIM-1; y++){
+    //             target[IX(0,x,y)] = source[IX(DIM-2,x,y)];
+    //         }
+    //     }
+    // } else {
         for(int x=1; x < DIM-1; x++){
             for(int y = 1; y < DIM-1; y++){
                 target[IX(0,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(1,x,y)] : target[IX(1,x,y)];
             }
         }
+    // }
+
+    // if(ARR_EXISTS(chunk_mask,2,1,1)){
+    //     source = GET_ARR(env,neighborArray,CK(2,1,1));
+    //     for(int x=1; x < DIM-1; x++){
+    //         for(int y = 1; y < DIM-1; y++){
+    //             target[IX(DIM-1,x,y)] = source[IX(1,x,y)];
+    //         }
+    //     }
+    // } else {
         for(int x=1; x < DIM-1; x++){
             for(int y = 1; y < DIM-1; y++){
                 target[IX(DIM-1,x,y)] =     vector_dir==BOUND_DIR_U ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)];
             }
         }
+    // }
     for(int x=1; x < DIM-1; x++){
         for(int y = 1; y < DIM-1; y++){
             //((x)+(DIM)*(y) + (DIM)*(DIM)*(z))
@@ -725,228 +598,3 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setBoundsToNeighbors
     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);
 }
-
-/**
- * This exclusively copies neighbors to make sure zeroing out stuff doesn't break sim
-*/
-JNIEXPORT void JNICALL Java_electrosphere_FluidSim_copyNeighbors
-  (JNIEnv * env,
-  jobject this,
-  jint N,
-  jint chunk_mask,
-  jint cx,
-  jint vector_dir,
-  jobjectArray neighborArray){
-    int DIM = N;
-    float * target = GET_ARR(env,neighborArray,CENTER_LOC);
-    float * source;
-
-
-    //
-    //
-    //   PLANES
-    //
-    //
-    if(ARR_EXISTS(chunk_mask,0,1,1)){
-        source = GET_ARR(env,neighborArray,CK(0,1,1));
-        for(int x=1; x < DIM-1; x++){
-            for(int y = 1; y < DIM-1; y++){
-                target[IX(0,x,y)] = source[IX(DIM-2,x,y)];
-            }
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,2,1,1)){
-        source = GET_ARR(env,neighborArray,CK(2,1,1));
-        for(int x=1; x < DIM-1; x++){
-            for(int y = 1; y < DIM-1; y++){
-                target[IX(DIM-1,x,y)] = source[IX(1,x,y)];
-            }
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,1,0,1)){
-        source = GET_ARR(env,neighborArray,CK(1,0,1));
-        for(int x=1; x < DIM-1; x++){
-            for(int y = 1; y < DIM-1; y++){
-                target[IX(x,0,y)] = source[IX(x,DIM-2,y)];
-            }
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,1,2,1)){
-        source = GET_ARR(env,neighborArray,CK(1,2,1));
-        for(int x=1; x < DIM-1; x++){
-            for(int y = 1; y < DIM-1; y++){
-                target[IX(x,DIM-1,y)] = source[IX(x,1,y)];
-            }
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,1,1,0)){
-        source = GET_ARR(env,neighborArray,CK(1,1,0));
-        for(int x=1; x < DIM-1; x++){
-            for(int y = 1; y < DIM-1; y++){
-                target[IX(x,y,0)] = source[IX(x,y,DIM-2)];
-            }
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,1,1,2)){
-        source = GET_ARR(env,neighborArray,CK(1,1,2));
-        for(int x=1; x < DIM-1; x++){
-            for(int y = 1; y < DIM-1; y++){
-                target[IX(x,y,DIM-1)] = source[IX(x,y,1)];
-            }
-        }
-    }
-
-
-    //
-    //
-    //    EDGES
-    //
-    //
-    if(ARR_EXISTS(chunk_mask,0,0,1)){
-        source = GET_ARR(env,neighborArray,CK(0,0,1));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(0,0,x)] = source[IX(DIM-2,DIM-2,x)];
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,2,0,1)){
-        source = GET_ARR(env,neighborArray,CK(2,0,1));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(DIM-1,0,x)] = source[IX(1,DIM-2,x)];
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,0,2,1)){
-        source = GET_ARR(env,neighborArray,CK(0,2,1));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(0,DIM-1,x)] = source[IX(DIM-2,1,x)];
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,2,2,1)){
-        source = GET_ARR(env,neighborArray,CK(2,2,1));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(DIM-1,DIM-1,x)] = source[IX(1,1,x)];
-        }
-    }
-
-    //
-    //
-
-    if(ARR_EXISTS(chunk_mask,0,1,0)){
-        source = GET_ARR(env,neighborArray,CK(0,1,0));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(0,x,0)] = source[IX(DIM-2,x,DIM-2)];
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,2,1,0)){
-        source = GET_ARR(env,neighborArray,CK(2,1,0));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(DIM-1,x,0)] = source[IX(1,x,DIM-2)];
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,0,1,2)){
-        source = GET_ARR(env,neighborArray,CK(0,1,2));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(0,x,DIM-1)] = source[IX(DIM-2,x,1)];
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,2,1,2)){
-        source = GET_ARR(env,neighborArray,CK(2,1,2));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(DIM-1,x,DIM-1)] = source[IX(1,x,1)];
-        }
-    }
-
-    //
-    //
-
-    if(ARR_EXISTS(chunk_mask,1,0,0)){
-        source = GET_ARR(env,neighborArray,CK(1,0,0));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(x,0,0)] = source[IX(x,DIM-2,DIM-2)];
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,1,2,0)){
-        source = GET_ARR(env,neighborArray,CK(1,2,0));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(x,DIM-1,0)] = source[IX(x,1,DIM-2)];
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,1,0,2)){
-        source = GET_ARR(env,neighborArray,CK(1,0,2));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(x,0,DIM-1)] = source[IX(x,DIM-2,1)];
-        }
-    }
-
-    if(ARR_EXISTS(chunk_mask,1,2,2)){
-        source = GET_ARR(env,neighborArray,CK(1,2,2));
-        for(int x=1; x < DIM-1; x++){
-                target[IX(x,DIM-1,DIM-1)] = source[IX(x,1,1)];
-        }
-    }
-
-
-    //
-    //
-    //     CORNERS
-    //
-    //
-
-    if(ARR_EXISTS(chunk_mask,0,0,0)){
-        source = GET_ARR(env,neighborArray,CK(0,0,0));
-        target[IX(0,0,0)] = source[IX(DIM-2,DIM-2,DIM-2)];
-    }
-
-    if(ARR_EXISTS(chunk_mask,2,0,0)){
-        source = GET_ARR(env,neighborArray,CK(2,0,0));
-        target[IX(DIM-1,0,0)] = source[IX(1,DIM-2,DIM-2)];
-    }
-
-    if(ARR_EXISTS(chunk_mask,0,2,0)){
-        source = GET_ARR(env,neighborArray,CK(0,2,0));
-        target[IX(0,DIM-1,0)] = source[IX(DIM-2,1,DIM-2)];
-    }
-
-    if(ARR_EXISTS(chunk_mask,2,2,0)){
-        source = GET_ARR(env,neighborArray,CK(2,2,0));
-        target[IX(DIM-1,DIM-1,0)] = source[IX(1,1,DIM-2)];
-    }
-
-    //
-    //
-
-    if(ARR_EXISTS(chunk_mask,0,0,2)){
-        source = GET_ARR(env,neighborArray,CK(0,0,2));
-        target[IX(0,0,DIM-1)] = source[IX(DIM-2,DIM-2,1)];
-    }
-
-    if(ARR_EXISTS(chunk_mask,2,0,2)){
-        source = GET_ARR(env,neighborArray,CK(2,0,2));
-        target[IX(DIM-1,0,DIM-1)] = source[IX(1,DIM-2,1)];
-    }
-
-    if(ARR_EXISTS(chunk_mask,0,2,2)){
-        source = GET_ARR(env,neighborArray,CK(0,2,2));
-        target[IX(0,DIM-1,DIM-1)] = source[IX(DIM-2,1,1)];
-    }
-
-    if(ARR_EXISTS(chunk_mask,2,2,2)){
-        source = GET_ARR(env,neighborArray,CK(2,2,2));
-        target[IX(DIM-1,DIM-1,DIM-1)] = source[IX(1,1,1)];
-    }
-
-
-
-}

src/main/java/electrosphere/FluidSim.java

@@ -42,7 +42,7 @@ public class FluidSim {
     //Buffers that contain density for current frame
     ByteBuffer[] density = new ByteBuffer[27];
     //Buffers that contain new density to add to the simulation
-    ByteBuffer[] densityAddition = new ByteBuffer[27];
+    ByteBuffer densityAddition;
     //Buffers that contain u vector directions
     ByteBuffer[] uVector = new ByteBuffer[27];
     //Buffers that contain v vector directions
@@ -66,23 +66,23 @@ public class FluidSim {
     float[] wArrayView = new float[DIM * DIM * DIM];
     //these should be set to the 
     float[] u0ArrayView = new float[DIM * DIM * DIM];
-    public float[] v0ArrayView = new float[DIM * DIM * DIM];
+    float[] v0ArrayView = new float[DIM * DIM * DIM];
     float[] w0ArrayView = new float[DIM * DIM * DIM];
 
     int chunkMask = 0;
     
 
-    static final float DIFFUSION_CONSTANT = 0.0f;
+    static final float DIFFUSION_CONSTANT = 0.00001f;
-    static final float VISCOSITY_CONSTANT = 0.0f;
+    static final float VISCOSITY_CONSTANT = 0.00001f;
 
-    static final int LINEARSOLVERTIMES = 20;
+    static final int LINEARSOLVERTIMES = 10;
 
-    static final float GRAVITY = -100f;
+    static final float GRAVITY = -1000f;
     
     public void setup(Vector3i offset){
         //allocate buffers for this chunk
         density[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-        densityAddition[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        densityAddition = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
         uVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
         vVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
         wVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
@@ -91,7 +91,7 @@ public class FluidSim {
         wAdditionVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
         //order endian-ness
         density[13].order(ByteOrder.LITTLE_ENDIAN);
-        densityAddition[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);
@@ -110,19 +110,19 @@ public class FluidSim {
         for(int i = 0; i < DIM; i++){
             for(int j = 0; j < DIM; j++){
                 for(int k = 0; k < DIM; k++){
-                    if(offset.x == 0 && offset.y == 0 && offset.z == 0){
+                    if(offset.x == 1){
                         if(
-                            Math.abs(16 - i) < 5 &&
+                            Math.abs(5 - i) < 4 &&
-                            Math.abs(j) < 5 &&
+                            Math.abs(j) < 4 &&
-                            Math.abs(16 - k) < 5 &&
+                            Math.abs(5 - k) < 4 &&
                             i < 17 && i > 0 &&
                             j < 17 && j > 0 &&
                             k < 17 && k > 0
                         ){
                             xf.put(1);
-                            uf.put(50);
+                            uf.put(1);
-                            vf.put(0);
+                            vf.put(-1f);
-                            wf.put(0);
+                            wf.put(rand.nextFloat() * 0.1f);
                         } else {
                             xf.put(0);
                             uf.put(0);
@@ -131,17 +131,17 @@ public class FluidSim {
                         }
                     } else {
                         if(
-                            Math.abs(0 - i) < 5 &&
+                            Math.abs(8 - i) < 4 &&
-                            Math.abs(j) < 5 &&
+                            Math.abs(j) < 4 &&
-                            Math.abs(0 - k) < 5 &&
+                            Math.abs(16 - k) < 4 &&
                             i < 17 && i > 0 &&
                             j < 17 && j > 0 &&
                             k < 17 && k > 0
                         ){
-                            // xf.put(1);
+                            xf.put(1);
-                            // uf.put(50);
+                            uf.put(1);
-                            // vf.put(0);
+                            vf.put(-1f);
-                            // wf.put(rand.nextFloat() * 0.1f);
+                            wf.put(rand.nextFloat() * 0.1f);
                         } else {
                             xf.put(0);
                             uf.put(0);
@@ -178,23 +178,20 @@ public class FluidSim {
         addVectorSources(simArray, timestep);
         swapAllVectorFields(simArray, timestep);
         solveVectorDiffusion(simArray, timestep);
-        solveProjection(simArray, step, timestep);
+        solveProjection(simArray, timestep);
         swapAllVectorFields(simArray, timestep);
         advectVectorsAcrossBoundaries(simArray, timestep);
-        solveProjection(simArray, step, timestep);
+        solveProjection(simArray, timestep);
 
         //
         //Density stage
         addDensity(simArray, timestep);
-        swapAllDensityArrays(simArray, timestep);
         diffuseDensity(simArray, timestep);
-        swapAllDensityArrays(simArray, timestep);
         advectDensity(simArray, timestep);
         // mirrorNeighborDensities(simArray, timestep);
         
 
 
-
         //
         //Read out of buffers back into accessible arrays
         for(int x = 0; x < simArray.length; x++){
@@ -209,128 +206,25 @@ public class FluidSim {
         }
     }
 
-    private static double sumAllDensity(FluidSim[][][] simArray){
+    /**
-        double rVal = 0;
+     * Runs a frame of the fluid simulation
-        for(int x = 0; x < simArray.length; x++){
+     */
-            for(int y = 0; y < simArray[0].length; y++){
+    private void simulate(int step, float timestep){
-                for(int z = 0; z < simArray[0][0].length; z++){
+        simulate(
-                    rVal = rVal + simArray[x][y][z].sumDensity();
+            DIM,
-                }
+            0,
-            }
+            density,
-        }
+            densityAddition,
-        return rVal;
+            uVector,
-    }
+            vVector,
-
+            wVector,
-    private double sumDensity(){
+            uAdditionVector,
-        double rVal = 0;
+            vAdditionVector,
-        for(int x = 1; x < DIM - 1; x++){
+            wAdditionVector,
-            for(int y = 1; y < DIM - 1; y++){
+            DIFFUSION_CONSTANT,
-                for(int z = 1; z < DIM - 1; z++){
+            VISCOSITY_CONSTANT,
-                    rVal = rVal + densityArrayView[IX(x,y,z)];
+            timestep
-                }
+        );
-            }
-        }
-        return rVal;
-    }
-
-    private double sumU(){
-        double rVal = 0;
-        for(int x = 1; x < DIM - 1; x++){
-            for(int y = 1; y < DIM - 1; y++){
-                for(int z = 1; z < DIM - 1; z++){
-                    rVal = rVal + Math.abs(uArrayView[IX(x,y,z)]);
-                }
-            }
-        }
-        return rVal;
-    }
-
-    private static double sumAllU(FluidSim[][][] simArray){
-        double rVal = 0;
-        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].readDataIntoArrays();
-                    rVal = rVal + simArray[x][y][z].sumU();
-                }
-            }
-        }
-        return rVal;
-    }
-
-    private double sumU0(){
-        double rVal = 0;
-        for(int x = 1; x < DIM - 1; x++){
-            for(int y = 1; y < DIM - 1; y++){
-                for(int z = 1; z < DIM - 1; z++){
-                    rVal = rVal + Math.abs(u0ArrayView[IX(x,y,z)]);
-                }
-            }
-        }
-        return rVal;
-    }
-
-    private static double sumAllU0(FluidSim[][][] simArray){
-        double rVal = 0;
-        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].readDataIntoArrays();
-                    rVal = rVal + simArray[x][y][z].sumU0();
-                }
-            }
-        }
-        return rVal;
-    }
-
-    private double sumV(){
-        double rVal = 0;
-        for(int x = 1; x < DIM - 1; x++){
-            for(int y = 1; y < DIM - 1; y++){
-                for(int z = 1; z < DIM - 1; z++){
-                    rVal = rVal + Math.abs(vArrayView[IX(x,y,z)]);
-                }
-            }
-        }
-        return rVal;
-    }
-
-    private static double sumAllV(FluidSim[][][] simArray){
-        double rVal = 0;
-        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].readDataIntoArrays();
-                    rVal = rVal + simArray[x][y][z].sumV();
-                }
-            }
-        }
-        return rVal;
-    }
-
-    private double sumV0(){
-        double rVal = 0;
-        for(int x = 1; x < DIM - 1; x++){
-            for(int y = 1; y < DIM - 1; y++){
-                for(int z = 1; z < DIM - 1; z++){
-                    rVal = rVal + Math.abs(v0ArrayView[IX(x,y,z)]);
-                }
-            }
-        }
-        return rVal;
-    }
-
-    private static double sumAllV0(FluidSim[][][] simArray){
-        double rVal = 0;
-        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].readDataIntoArrays();
-                    rVal = rVal + simArray[x][y][z].sumV0();
-                }
-            }
-        }
-        return rVal;
     }
 
     private static void solveChunkMask(FluidSim[][][] simArray){
@@ -362,20 +256,6 @@ public class FluidSim {
     }
 
     private static void solveVectorDiffusion(FluidSim[][][] simArray, float timestep){
-        //samples u,v,w,u0,v0,w0
-        //sets u,v,w
-        // 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].copyNeighborsWrapper(1, simArray[x][y][z].uVector);
-        //             simArray[x][y][z].copyNeighborsWrapper(2, simArray[x][y][z].vVector);
-        //             simArray[x][y][z].copyNeighborsWrapper(3, simArray[x][y][z].wVector);
-        //             simArray[x][y][z].copyNeighborsWrapper(1, simArray[x][y][z].uAdditionVector);
-        //             simArray[x][y][z].copyNeighborsWrapper(2, simArray[x][y][z].vAdditionVector);
-        //             simArray[x][y][z].copyNeighborsWrapper(3, simArray[x][y][z].wAdditionVector);
-        //         }
-        //     }
-        // }
         for(int l = 0; l < LINEARSOLVERTIMES; l++){
             for(int x = 0; x < simArray.length; x++){
                 for(int y = 0; y < simArray[0].length; y++){
@@ -394,38 +274,16 @@ public class FluidSim {
                         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);
-                        simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
-                        simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
-                        simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
                     }
                 }
             }
         }
     }
 
-    private static void solveProjection(FluidSim[][][] simArray, int step, float timestep){
+    private static void solveProjection(FluidSim[][][] simArray, float timestep){
-        //samples u,v,w
-        //sets u0,v0
         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);
-                    simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uAdditionVector);
-                    simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
-                    simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
-                    simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
-                    simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vAdditionVector);
-                }
-            }
-        }
-        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++){
-                    // System.out.println("Setup " + x + " " +  y + " " + z);
                     //setup projection across boundaries
                     //...
                     //set boundaries appropriately
@@ -439,14 +297,9 @@ public class FluidSim {
                 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);
-                    simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].uAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].vAdditionVector);
                 }
             }
         }
-        //samples u0, v0
-        //sets u0
-        //these should have just been mirrored in the above
         //
         //Perform main projection solver
         for(int l = 0; l < LINEARSOLVERTIMES; l++){
@@ -464,13 +317,11 @@ public class FluidSim {
                 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].copyNeighborsWrapper(0, x, simArray[x][y][z].uAdditionVector);
+                        // simArray[x][y][z].setBoundsToNeighborsWrapper(0, simArray[x][y][z].vAdditionVector);
                     }
                 }
             }
         }
-        //samples u,v,w,u0
-        //sets u,v,w
         //Finalize projection
         for(int x = 0; x < simArray.length; x++){
             for(int y = 0; y < simArray[0].length; y++){
@@ -489,15 +340,6 @@ public class FluidSim {
                     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);
-                    simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
-                    simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
-                    simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
-                    simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uAdditionVector);
-                    simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vAdditionVector);
-                    simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].uAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].vAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].wAdditionVector);
                 }
             }
         }
@@ -521,19 +363,6 @@ public class FluidSim {
                 }
             }
         }
-        //then need to mirror each array as relevant
-        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].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
-                    simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
-                    simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
-                    simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].wAdditionVector);
-                }
-            }
-        }
     }
 
     private void swapVectorFields(){
@@ -543,42 +372,21 @@ public class FluidSim {
         // densityAddition = density[13];
         // density[13] = tmp;
         //swap u0 <-> u
-        for(int i = 0; i < 27; i++){
+        tmp = uAdditionVector[13];
-            tmp = uAdditionVector[i];
+        uAdditionVector[13] = uVector[13];
-            uAdditionVector[i] = uVector[i];
+        uVector[13] = tmp;
-            uVector[i] = tmp;
         //swap v0 <-> v
-            tmp = vAdditionVector[i];
+        tmp = vAdditionVector[13];
-            vAdditionVector[i] = vVector[i];
+        vAdditionVector[13] = vVector[13];
-            vVector[i] = tmp;
+        vVector[13] = tmp;
         //swap w0 <-> w
-            tmp = wAdditionVector[i];
+        tmp = wAdditionVector[13];
-            wAdditionVector[i] = wVector[i];
+        wAdditionVector[13] = wVector[13];
-            wVector[i] = tmp;
+        wVector[13] = tmp;
-        }
         //...
     }
 
     private static void advectVectorsAcrossBoundaries(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++){
-                    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);
-                    simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uAdditionVector);
-                    simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vAdditionVector);
-                    simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
-                    simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
-                    simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
-                    simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vAdditionVector);
-                    simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wAdditionVector);
-                }
-            }
-        }
-        //samples u,v,w,u0,v0,w0
         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++){
@@ -597,9 +405,6 @@ public class FluidSim {
                     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);
-                    simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
-                    simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
-                    simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
                 }
             }
         }
@@ -613,36 +418,16 @@ public class FluidSim {
                     simArray[x][y][z].addDensityWrapper(timestep);
                     //swap x <=> x0
                     //swap arrays in java side...
-                    // simArray[x][y][z].swapDensityArrays();
-                }
-            }
-        }
-    }
-
-    private static void swapAllDensityArrays(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++){
                     simArray[x][y][z].swapDensityArrays();
                 }
             }
         }
-        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].copyNeighborsWrapper(0, x, simArray[x][y][z].density);
-                    simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].densityAddition);
-                }
-            }
-        }
     }
 
     private void swapDensityArrays(){
-        for(int i = 0; i < 27; i++){
+        ByteBuffer tmp = density[13];
-            ByteBuffer tmp = density[i];
+        density[13] = densityAddition;
-            density[i] = densityAddition[i];
+        densityAddition = tmp;
-            densityAddition[i] = tmp;
-        }
     }
 
     private static void diffuseDensity(FluidSim[][][] simArray, float timestep){
@@ -672,7 +457,7 @@ public class FluidSim {
             for(int y = 0; y < simArray[0].length; y++){
                 for(int z = 0; z < simArray[0][0].length; z++){
                     //swap x <=> x0 again
-                    // simArray[x][y][z].swapDensityArrays();
+                    simArray[x][y][z].swapDensityArrays();
                     //advect density
                     simArray[x][y][z].advectDensityWrapper(timestep);
                 }
@@ -790,6 +575,14 @@ public class FluidSim {
     }
     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
      */
@@ -812,7 +605,7 @@ public class FluidSim {
     private void addDensityWrapper(float timestep){
         addDensity(DIM, chunkMask, density, densityAddition, timestep);
     }
-    private native void addDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer[] x0, float timestep);
+    private native void addDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer x0, float timestep);
 
     /**
      * Solve density diffusion
@@ -820,7 +613,7 @@ public class FluidSim {
     private void solveDiffuseDensityWrapper(float timestep){
         solveDiffuseDensity(DIM, chunkMask, density, densityAddition, uVector, vVector, wVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
     }
-    private native void solveDiffuseDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer[] x0, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+    private native void solveDiffuseDensity(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
@@ -828,7 +621,7 @@ public class FluidSim {
     private void advectDensityWrapper(float timestep){
         advectDensity(DIM, chunkMask, density, densityAddition, uVector, vVector, wVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
     }
-    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);
+    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);
 
 
     /**
@@ -839,16 +632,6 @@ public class FluidSim {
     }
     private native void setBoundsToNeighbors(int DIM_X, int chunkMask, int vectorDir, ByteBuffer[] neighborMap);
 
-    /**
-     * Sets the bounds of the neighbormap to neighbor values if available, otherwise doesn't mess with them.
-     * This is to make sure zeroing out doesn't mess up the sim
-     */
-    private void copyNeighborsWrapper(int vectorDir, int x, ByteBuffer[] neighborMap){
-        copyNeighbors(DIM, chunkMask, x, vectorDir, neighborMap);
-    }
-    private native void copyNeighbors(int DIM_X, int chunkMask, int x, int vectorDir, ByteBuffer[] neighborMap);
-
-
 
 
 
@@ -920,34 +703,17 @@ public class FluidSim {
         if(wVector[13].position() > 0){
             wVector[13].position(0);
         }
-        if(uAdditionVector[13].position() > 0){
-            uAdditionVector[13].position(0);
-        }
-        if(vAdditionVector[13].position() > 0){
-            vAdditionVector[13].position(0);
-        }
-        if(wAdditionVector[13].position() > 0){
-            wAdditionVector[13].position(0);
-        }
         FloatBuffer xFloatView = density[13].asFloatBuffer();
         FloatBuffer uFloatView = uVector[13].asFloatBuffer();
         FloatBuffer vFloatView = vVector[13].asFloatBuffer();
         FloatBuffer wFloatView = wVector[13].asFloatBuffer();
-        FloatBuffer u0FloatView = uAdditionVector[13].asFloatBuffer();
-        FloatBuffer v0FloatView = vAdditionVector[13].asFloatBuffer();
-        FloatBuffer w0FloatView = wAdditionVector[13].asFloatBuffer();
-        int index = 0;
         for(int i = 0; i < DIM; i++){
             for(int j = 0; j < DIM; j++){
                 for(int k = 0; k < DIM; k++){
-                    index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
+                    densityArrayView[IX(i,j,k)] = xFloatView.get();
-                    densityArrayView[index] = xFloatView.get();
+                    uArrayView[IX(i,j,k)] = uFloatView.get();
-                    uArrayView[index] = uFloatView.get();
+                    vArrayView[IX(i,j,k)] = vFloatView.get();
-                    vArrayView[index] = vFloatView.get();
+                    wArrayView[IX(i,j,k)] = wFloatView.get();
-                    wArrayView[index] = wFloatView.get();
-                    u0ArrayView[index] = u0FloatView.get();
-                    v0ArrayView[index] = v0FloatView.get();
-                    w0ArrayView[index] = w0FloatView.get();
                 }
             }
         }
@@ -957,8 +723,8 @@ public class FluidSim {
      * Writes data from the java-side arrays into buffers that get passed into c-side
      */
     private void writeNewStateIntoBuffers(){
-        if(densityAddition[13].position() > 0){
+        if(densityAddition.position() > 0){
-            densityAddition[13].position(0);
+            densityAddition.position(0);
         }
         if(uAdditionVector[13].position() > 0){
             uAdditionVector[13].position(0);
@@ -969,19 +735,17 @@ public class FluidSim {
         if(wAdditionVector[13].position() > 0){
             wAdditionVector[13].position(0);
         }
-        FloatBuffer x0FloatView = densityAddition[13].asFloatBuffer();
+        FloatBuffer x0FloatView = densityAddition.asFloatBuffer();
         FloatBuffer u0FloatView = uAdditionVector[13].asFloatBuffer();
         FloatBuffer v0FloatView = vAdditionVector[13].asFloatBuffer();
         FloatBuffer w0FloatView = wAdditionVector[13].asFloatBuffer();
-        int index = 0;
         for(int i = 0; i < DIM; i++){
             for(int j = 0; j < DIM; j++){
                 for(int k = 0; k < DIM; k++){
-                    index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
+                    x0FloatView.put(density0ArrayView[IX(i,j,k)]);
-                    x0FloatView.put(density0ArrayView[index]);
+                    u0FloatView.put(u0ArrayView[IX(i,j,k)]);
-                    u0FloatView.put(u0ArrayView[index]);
+                    v0FloatView.put(v0ArrayView[IX(i,j,k)]);
-                    v0FloatView.put(v0ArrayView[index]);
+                    w0FloatView.put(w0ArrayView[IX(i,j,k)]);
-                    w0FloatView.put(w0ArrayView[index]);
                 }
             }
         }
@@ -991,14 +755,10 @@ public class FluidSim {
      * Adds gravity to the simulation
      */
     private void addGravity(){
-        int index = 0;
         for(int i = 0; i < DIM; i++){
             for(int j = 0; j < DIM; j++){
                 for(int k = 0; k < DIM; k++){
-                    index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
+                    v0ArrayView[IX(i,j,k)] = densityArrayView[IX(i,j,k)] * GRAVITY;
-                    u0ArrayView[index] = 0;
-                    v0ArrayView[index] = densityArrayView[index] * GRAVITY;
-                    w0ArrayView[index] = 0;
                 }
             }
         }
@@ -1063,10 +823,6 @@ public class FluidSim {
         return density[getNeighborIndex(1,1,1)];
     }
 
-    public ByteBuffer getDensityAdditionBuffer(){
-        return densityAddition[getNeighborIndex(1,1,1)];
-    }
-
     public ByteBuffer getUBuffer(){
         return uVector[getNeighborIndex(1,1,1)];
     }
@@ -1093,7 +849,6 @@ public class FluidSim {
 
     public void setNeighbor(int x, int y, int z, FluidSim neighbor){
         density[getNeighborIndex(x,y,z)] = neighbor.getDensityBuffer();
-        densityAddition[getNeighborIndex(x,y,z)] = neighbor.getDensityAdditionBuffer();
         uVector[getNeighborIndex(x,y,z)] = neighbor.getUBuffer();
         vVector[getNeighborIndex(x,y,z)] = neighbor.getVBuffer();
         wVector[getNeighborIndex(x,y,z)] = neighbor.getWBuffer();

src/main/java/electrosphere/Main.java

@@ -22,7 +22,7 @@ public class Main {
 
     public static void main(String args[]){
 
-        int dim = 10;
+        int dim = 2;
         int i = 0;
         long time = 0;
         long lastTime = 0;
@@ -55,7 +55,7 @@ public class Main {
                 //
                 //Simulate
                 //
-                FluidSim.simChunks(simArray,i,0.01f);
+                FluidSim.simChunks(simArray,i,0.001f);
                 time = time + (System.currentTimeMillis() - lastTime);
                 //
                 //Remesh
@@ -70,8 +70,8 @@ public class Main {
                 //redraw
                 GLFWContext.redraw(meshArray);
                 i++;
-                if(i == 100){
+                if(i == 1000){
-                    System.out.println(time / 100.0);
+                    System.out.println(time / 1000.0);
                 }
                 if(i > 3){
                     // scan.next();

src/main/java/electrosphere/render/Mesh.java

@@ -2,13 +2,11 @@ package electrosphere.render;
 
 import java.io.BufferedInputStream;
 import java.io.BufferedReader;
-import java.io.File;
 import java.io.IOException;
 import java.io.InputStream;
 import java.io.InputStreamReader;
 import java.nio.FloatBuffer;
 import java.nio.IntBuffer;
-import java.nio.file.Files;
 import java.util.HashMap;
 import java.util.LinkedList;
 import java.util.List;
@@ -243,7 +241,7 @@ public class Mesh {
     public static void initShaderProgram(){
         String vsSrc = "";
         ClassLoader classloader = Thread.currentThread().getContextClassLoader();
-        try (BufferedReader is = new BufferedReader(Files.newBufferedReader(new File("C:\\Users\\satellite\\Documents\\fluid-sim\\src\\main\\resources\\shader.vs").toPath()))){
+        try (BufferedReader is = new BufferedReader(new InputStreamReader(classloader.getResourceAsStream("shader.vs")))){
             String temp;
             while((temp = is.readLine())!=null){
                 vsSrc = vsSrc + temp + "\n";
@@ -254,7 +252,7 @@ public class Mesh {
         }
 
         String fsSrc = "";
-        try (BufferedReader is = new BufferedReader(Files.newBufferedReader(new File("C:\\Users\\satellite\\Documents\\fluid-sim\\src\\main\\resources\\shader.fs").toPath()))){
+        try (BufferedReader is = new BufferedReader(new InputStreamReader(classloader.getResourceAsStream("shader.fs")))){
             String temp;
             while((temp = is.readLine())!=null){
                 fsSrc = fsSrc + temp + "\n";