update jenkinsfile

.gitignore

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

.gitmodules

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

Jenkinsfile

@@ -0,0 +1,17 @@
+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/compile.sh

@@ -42,22 +42,17 @@ rm -f ./*.dll
 
 #compile object files
 COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
-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"
+INPUT_FILES="./src/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="./velocitystep.c"
+INPUT_FILES="./src/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="./chunkmask.c"
+INPUT_FILES="./src/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/fluidsim.c

@@ -1,396 +0,0 @@
-#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 1.0E-5f
+#define electrosphere_FluidSim_DIFFUSION_CONSTANT 0.0f
 #undef electrosphere_FluidSim_VISCOSITY_CONSTANT
-#define electrosphere_FluidSim_VISCOSITY_CONSTANT 1.0E-5f
+#define electrosphere_FluidSim_VISCOSITY_CONSTANT 0.0f
 #undef electrosphere_FluidSim_LINEARSOLVERTIMES
-#define electrosphere_FluidSim_LINEARSOLVERTIMES 10L
+#define electrosphere_FluidSim_LINEARSOLVERTIMES 20L
 #undef electrosphere_FluidSim_GRAVITY
-#define electrosphere_FluidSim_GRAVITY -1000.0f
+#define electrosphere_FluidSim_GRAVITY -100.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, jobject, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, 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, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, 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, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
@@ -113,6 +113,14 @@ 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

@@ -0,0 +1,10 @@
+//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

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

src/main/c/src/chunkmask.c

@@ -1,7 +1,8 @@
 #include <jni.h>
 #include <stdint.h>
-#include "includes/utilities.h"
-#include "includes/chunkmask.h"
+#include "../includes/libfluidsim.h"
+#include "../includes/utilities.h"
+#include "../includes/chunkmask.h"
 
 uint32_t matrix_transform(JNIEnv * env, jobjectArray jrx);
 
@@ -60,43 +61,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/src/densitystep.c

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

src/main/c/src/velocitystep.c

@@ -3,8 +3,9 @@
 #include <immintrin.h>
 #include <stdint.h>
 
-#include "includes/utilities.h"
-#include "includes/chunkmask.h"
+#include "../includes/libfluidsim.h"
+#include "../includes/utilities.h"
+#include "../includes/chunkmask.h"
 
 
 #define BOUND_NO_DIR 0
@@ -16,7 +17,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, float * d0, float * u, float * v, float * w, 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);
 
 
 /*
@@ -66,7 +67,7 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveVectorDiffuse
   jfloat DIFFUSION_CONST,
   jfloat VISCOSITY_CONST,
   jfloat dt){
-    float a=dt*DIFFUSION_CONST*N*N*N;
+    float a=dt*VISCOSITY_CONST*N*N*N;
     float c=1+6*a;
     int i, j, k, l, m;
     float * u = GET_ARR(env,jru,CENTER_LOC);
@@ -177,7 +178,9 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupProjection
   jfloat dt){
     int i, j, k;
 
-    __m256 nVector = _mm256_set1_ps(N);
+    __m256 xVector = _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;
@@ -189,6 +192,9 @@ 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;
@@ -198,15 +204,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,nVector);
+            vector = _mm256_div_ps(vector,xVector);
             //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);
+            vector2 = _mm256_div_ps(vector2,yVector);
             //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);
+            vector3 = _mm256_div_ps(vector3,zVector);
             //multiply and finalize
             vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
             vector = _mm256_mul_ps(vector,constScalar);
@@ -220,21 +226,30 @@ 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,nVector);
+            vector = _mm256_div_ps(vector,xVector);
             //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);
+            vector2 = _mm256_div_ps(vector2,yVector);
             //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);
+            vector3 = _mm256_div_ps(vector3,zVector);
             //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;
+            // }
         }
     }
 }
@@ -262,31 +277,34 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveProjection
     __m256 aScalar = _mm256_set1_ps(a);
     __m256 cScalar = _mm256_set1_ps(c);
 
-    float * x = GET_ARR(env,jru0,CENTER_LOC);
-    float * x0 = GET_ARR(env,jrv0,CENTER_LOC);
+    float * p = GET_ARR(env,jru0,CENTER_LOC);
+    float * div = GET_ARR(env,jrv0,CENTER_LOC);
     // update for each cell
     for(k=1; k<N-1; k++){
         for(j=1; j<N-1; j++){
             int n = 0;
             //solve as much as possible vectorized
             for(i = 1; i < N-1; i=i+8){
-                __m256 vector = _mm256_loadu_ps(&x[IX(i-1,j,k)]);
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i+1,j,k)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j-1,k)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j+1,k)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k-1)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k+1)]));
-                vector = _mm256_mul_ps(vector,aScalar);
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x0[IX(i,j,k)]));
+                __m256 vector = _mm256_loadu_ps(&p[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(&p[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(&p[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j,k+1)]));
+                // vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&div[IX(i,j,k)]));
                 vector = _mm256_div_ps(vector,cScalar);
-                _mm256_storeu_ps(&x[IX(i,j,k)],vector);
+                _mm256_storeu_ps(&p[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;
+                    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;
                 }
             }
+            // 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;
+            // }
         }
     }
 }
@@ -309,9 +327,10 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection
   jfloat VISCOSITY_CONST,
   jfloat dt){
     int i, j, k;
-    __m256 nVector = _mm256_set1_ps(N);
-    __m256 constScalar = _mm256_set1_ps(0.5f*N);
-    __m256 zeroVec = _mm256_set1_ps(0);
+    // __m256 constScalar = _mm256_set1_ps(0.5f*N);
+    __m256 xScalar = _mm256_set1_ps(0.5*N);
+    __m256 yScalar = _mm256_set1_ps(0.5*N);
+    __m256 zScalar = _mm256_set1_ps(0.5*N);
     __m256 vector, vector2, vector3;
 
     float * u = GET_ARR(env,jru,CENTER_LOC);
@@ -321,6 +340,8 @@ 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++ ) {
             //
@@ -330,14 +351,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,constScalar);
+            vector = _mm256_mul_ps(vector,xScalar);
             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_mul_ps(vector,xScalar);
             vector = _mm256_sub_ps(_mm256_loadu_ps(&u[IX(9,j,k)]),vector);
             _mm256_storeu_ps(&u[IX(9,j,k)],vector);
             //
@@ -347,14 +368,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,constScalar);
+            vector = _mm256_mul_ps(vector,yScalar);
             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_mul_ps(vector,yScalar);
             vector = _mm256_sub_ps(_mm256_loadu_ps(&v[IX(9,j,k)]),vector);
             _mm256_storeu_ps(&v[IX(9,j,k)],vector);
             //
@@ -364,16 +385,21 @@ 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,constScalar);
+            vector = _mm256_mul_ps(vector,zScalar);
             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_mul_ps(vector,zScalar);
             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;
+            // }
         }
     }
 }
@@ -395,13 +421,13 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectVectors
   jfloat DIFFUSION_CONST,
   jfloat VISCOSITY_CONST,
   jfloat dt){
-    advect(env,chunk_mask,N,1,jru,GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
-    advect(env,chunk_mask,N,2,jrv,GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
-    advect(env,chunk_mask,N,3,jrw,GET_ARR(env,jrw0,CENTER_LOC),GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
+    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,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,3,jrw,jrw0,GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
 }
 
 
-void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt){
+void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, jobjectArray jrd0, 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;
@@ -410,12 +436,12 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, f
 
     float * d = GET_ARR(env,jrd,CENTER_LOC);
 
-    float * sampleArr = d0;
+    float * d0 = GET_ARR(env,jrd0,CENTER_LOC);
 
     for(k=1; k<N-1; k++){
         for(j=1; j<N-1; j++){
             for(i=1; i<N-1; i++){
-                sampleArr = d0;
+                d0 = GET_ARR(env,jrd0,CENTER_LOC);
                 //calculate location to pull from
                 x = i-dtx*u[IX(i,j,k)];
                 y = j-dty*v[IX(i,j,k)];
@@ -423,55 +449,175 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, f
 
                 m = n = o = 1;
 
-                if(x < 1){ m -= 1; }
-                if(x >= N-1){ m += 1; }
-                if(y < 1){ n -= 1; }
-                if(y >= N-1){ n += 1; }
-                if(z < 1){ o -= 1; }
-                if(z >= N-1){ o += 1; }
+                if(x < 0){ m += 1; }
+                else if(x >= N){ m -= 1; }
+                if(y < 0){ n += 1; }
+                else if(y >= N){ n -= 1; }
+                if(z < 0){ o += 1; }
+                else if(z >= N){ o -= 1; }
 
                 //If the out of bounds coordinate is in bounds for a neighbor chunk, use that chunk as source instead
-                // if(CK(m,n,o) != CENTER_LOC){
-                //     printf("Looking in border chunk\n");
-                // }
-                // if(x > 16){
-                //     printf("%f %d %d %d\n",m,n,o);
-                // }
-                // if(CK(m,n,o) != CENTER_LOC && ARR_EXISTS(chunk_mask,m,n,o)){
-                //     // printf("Hit other chunk\n");
-                //     sampleArr = GET_ARR(env,jrd,CK(m,n,o));
-                //     x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * N;
-                //     y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * N;
-                //     z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * N;
+                if(CK(m,n,o) != CENTER_LOC && ARR_EXISTS(chunk_mask,m,n,o)){
+
+                    // if(i == 1 && j == 1 && k == 1){
+                    //     printf("\narr indices: %d %d %d\n\n",m,n,o);
                     // }
 
+                    //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
-                if (x<0.5f) x=0.5f;
-                if (x>N+0.5f) x=N+0.5f;
-                if (y<0.5f) y=0.5f;
-                if (y>N+0.5f) y=N+0.5f;
-                if (z<0.5f) z=0.5f;
-                if (z>N+0.5f) z=N+0.5f;
-
-                //get actual indices
+                //get indices, and calculate percentage to pull from each index
+                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;
-                j0=(int)y;
-                j1=j0+1;
-                k0=(int)z;
-                k1=k0+1;
-
-                //calculate percentage of each index
                     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(i0 >= N){
-                    i0 = N - 1;
                 }
+
+                // 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;
+                // i1=i0+1;
+                // j0=(int)y;
+                // j1=j0+1;
+                // k0=(int)z;
+                // k1=k0+1;
+
+                //calculate percentage of each index
+                // s1 = x-i0;
+                // s0 = 1-s1;
+                // t1 = y-j0;
+                // t0 = 1-t1;
+                // u1 = z-k0;
+                // u0 = 1-u1;
+
+                // if(i0 >= N){
+                //     i0 = N - 1;
+                // }
                 // if(i0 < 0){
                 //     i0 = 0;
                 // }
@@ -487,9 +633,9 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, f
                 // if(k0 < 0){
                 //     k0 = 0;
                 // }
-                if(i1 >= N){
-                    i1 = N - 1;
-                }
+                // if(i1 >= N){
+                //     i1 = N - 1;
+                // }
                 // if(i1 < 0){
                 //     i1 = 0;
                 // }
@@ -507,16 +653,16 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, f
                 // }
                 d[IX(i,j,k)] = 
                 s0*(
-                    t0*u0*sampleArr[IX(i0,j0,k0)]+
-                    t1*u0*sampleArr[IX(i0,j1,k0)]+
-                    t0*u1*sampleArr[IX(i0,j0,k1)]+
-                    t1*u1*sampleArr[IX(i0,j1,k1)]
+                    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*sampleArr[IX(i1,j0,k0)]+
-                    t1*u0*sampleArr[IX(i1,j1,k0)]+
-                    t0*u1*sampleArr[IX(i1,j0,k1)]+
-                    t1*u1*sampleArr[IX(i1,j1,k1)]
+                    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)]
                 );
             }
         }
@@ -533,35 +679,16 @@ 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))
@@ -598,3 +725,228 @@ 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;
+    ByteBuffer[] densityAddition = new ByteBuffer[27];
     //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];
-    float[] v0ArrayView = new float[DIM * DIM * DIM];
+    public float[] v0ArrayView = new float[DIM * DIM * DIM];
     float[] w0ArrayView = new float[DIM * DIM * DIM];
 
     int chunkMask = 0;
     
 
-    static final float DIFFUSION_CONSTANT = 0.00001f;
-    static final float VISCOSITY_CONSTANT = 0.00001f;
+    static final float DIFFUSION_CONSTANT = 0.0f;
+    static final float VISCOSITY_CONSTANT = 0.0f;
 
-    static final int LINEARSOLVERTIMES = 10;
+    static final int LINEARSOLVERTIMES = 20;
 
-    static final float GRAVITY = -1000f;
+    static final float GRAVITY = -100f;
     
     public void setup(Vector3i offset){
         //allocate buffers for this chunk
         density[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-        densityAddition = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        densityAddition[13] = 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.order(ByteOrder.LITTLE_ENDIAN);
+        densityAddition[13].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 == 1){
+                    if(offset.x == 0 && offset.y == 0 && offset.z == 0){
                         if(
-                            Math.abs(5 - i) < 4 &&
-                            Math.abs(j) < 4 &&
-                            Math.abs(5 - k) < 4 &&
+                            Math.abs(16 - i) < 5 &&
+                            Math.abs(j) < 5 &&
+                            Math.abs(16 - k) < 5 &&
                             i < 17 && i > 0 &&
                             j < 17 && j > 0 &&
                             k < 17 && k > 0
                         ){
                             xf.put(1);
-                            uf.put(1);
-                            vf.put(-1f);
-                            wf.put(rand.nextFloat() * 0.1f);
+                            uf.put(50);
+                            vf.put(0);
+                            wf.put(0);
                         } else {
                             xf.put(0);
                             uf.put(0);
@@ -131,17 +131,17 @@ public class FluidSim {
                         }
                     } else {
                         if(
-                            Math.abs(8 - i) < 4 &&
-                            Math.abs(j) < 4 &&
-                            Math.abs(16 - k) < 4 &&
+                            Math.abs(0 - i) < 5 &&
+                            Math.abs(j) < 5 &&
+                            Math.abs(0 - k) < 5 &&
                             i < 17 && i > 0 &&
                             j < 17 && j > 0 &&
                             k < 17 && k > 0
                         ){
-                            xf.put(1);
-                            uf.put(1);
-                            vf.put(-1f);
-                            wf.put(rand.nextFloat() * 0.1f);
+                            // xf.put(1);
+                            // uf.put(50);
+                            // vf.put(0);
+                            // wf.put(rand.nextFloat() * 0.1f);
                         } else {
                             xf.put(0);
                             uf.put(0);
@@ -178,20 +178,23 @@ public class FluidSim {
         addVectorSources(simArray, timestep);
         swapAllVectorFields(simArray, timestep);
         solveVectorDiffusion(simArray, timestep);
-        solveProjection(simArray, timestep);
+        solveProjection(simArray, step, timestep);
         swapAllVectorFields(simArray, timestep);
         advectVectorsAcrossBoundaries(simArray, timestep);
-        solveProjection(simArray, timestep);
+        solveProjection(simArray, step, 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++){
@@ -206,25 +209,128 @@ public class FluidSim {
         }
     }
 
-    /**
-     * Runs a frame of the fluid simulation
-     */
-    private void simulate(int step, float timestep){
-        simulate(
-            DIM,
-            0,
-            density,
-            densityAddition,
-            uVector,
-            vVector,
-            wVector,
-            uAdditionVector,
-            vAdditionVector,
-            wAdditionVector,
-            DIFFUSION_CONSTANT,
-            VISCOSITY_CONSTANT,
-            timestep
-        );
+    private static double sumAllDensity(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++){
+                    rVal = rVal + simArray[x][y][z].sumDensity();
+                }
+            }
+        }
+        return rVal;
+    }
+
+    private double sumDensity(){
+        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 + densityArrayView[IX(x,y,z)];
+                }
+            }
+        }
+        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){
@@ -256,6 +362,20 @@ 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++){
@@ -274,16 +394,38 @@ 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, float timestep){
+    private static void solveProjection(FluidSim[][][] simArray, int step, 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
@@ -297,9 +439,14 @@ 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++){
@@ -317,11 +464,13 @@ 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].setBoundsToNeighborsWrapper(0, simArray[x][y][z].vAdditionVector);
+                        simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].uAdditionVector);
                     }
                 }
             }
         }
+        //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++){
@@ -340,6 +489,15 @@ 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);
                 }
             }
         }
@@ -363,6 +521,19 @@ 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(){
@@ -372,21 +543,42 @@ public class FluidSim {
         // densityAddition = density[13];
         // density[13] = tmp;
         //swap u0 <-> u
-        tmp = uAdditionVector[13];
-        uAdditionVector[13] = uVector[13];
-        uVector[13] = tmp;
+        for(int i = 0; i < 27; i++){
+            tmp = uAdditionVector[i];
+            uAdditionVector[i] = uVector[i];
+            uVector[i] = tmp;
             //swap v0 <-> v
-        tmp = vAdditionVector[13];
-        vAdditionVector[13] = vVector[13];
-        vVector[13] = tmp;
+            tmp = vAdditionVector[i];
+            vAdditionVector[i] = vVector[i];
+            vVector[i] = tmp;
             //swap w0 <-> w
-        tmp = wAdditionVector[13];
-        wAdditionVector[13] = wVector[13];
-        wVector[13] = tmp;
+            tmp = wAdditionVector[i];
+            wAdditionVector[i] = wVector[i];
+            wVector[i] = 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++){
@@ -405,6 +597,9 @@ 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);
                 }
             }
         }
@@ -418,16 +613,36 @@ 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(){
-        ByteBuffer tmp = density[13];
-        density[13] = densityAddition;
-        densityAddition = tmp;
+        for(int i = 0; i < 27; i++){
+            ByteBuffer tmp = density[i];
+            density[i] = densityAddition[i];
+            densityAddition[i] = tmp;
+        }
     }
 
     private static void diffuseDensity(FluidSim[][][] simArray, float timestep){
@@ -457,7 +672,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);
                 }
@@ -575,14 +790,6 @@ 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
      */
@@ -605,7 +812,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
@@ -613,7 +820,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
@@ -621,7 +828,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);
 
 
     /**
@@ -632,6 +839,16 @@ 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);
+
+
 
 
 
@@ -703,17 +920,34 @@ 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++){
-                    densityArrayView[IX(i,j,k)] = xFloatView.get();
-                    uArrayView[IX(i,j,k)] = uFloatView.get();
-                    vArrayView[IX(i,j,k)] = vFloatView.get();
-                    wArrayView[IX(i,j,k)] = wFloatView.get();
+                    index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
+                    densityArrayView[index] = xFloatView.get();
+                    uArrayView[index] = uFloatView.get();
+                    vArrayView[index] = vFloatView.get();
+                    wArrayView[index] = wFloatView.get();
+                    u0ArrayView[index] = u0FloatView.get();
+                    v0ArrayView[index] = v0FloatView.get();
+                    w0ArrayView[index] = w0FloatView.get();
                 }
             }
         }
@@ -723,8 +957,8 @@ public class FluidSim {
      * Writes data from the java-side arrays into buffers that get passed into c-side
      */
     private void writeNewStateIntoBuffers(){
-        if(densityAddition.position() > 0){
-            densityAddition.position(0);
+        if(densityAddition[13].position() > 0){
+            densityAddition[13].position(0);
         }
         if(uAdditionVector[13].position() > 0){
             uAdditionVector[13].position(0);
@@ -735,17 +969,19 @@ public class FluidSim {
         if(wAdditionVector[13].position() > 0){
             wAdditionVector[13].position(0);
         }
-        FloatBuffer x0FloatView = densityAddition.asFloatBuffer();
+        FloatBuffer x0FloatView = densityAddition[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++){
-                    x0FloatView.put(density0ArrayView[IX(i,j,k)]);
-                    u0FloatView.put(u0ArrayView[IX(i,j,k)]);
-                    v0FloatView.put(v0ArrayView[IX(i,j,k)]);
-                    w0FloatView.put(w0ArrayView[IX(i,j,k)]);
+                    index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
+                    x0FloatView.put(density0ArrayView[index]);
+                    u0FloatView.put(u0ArrayView[index]);
+                    v0FloatView.put(v0ArrayView[index]);
+                    w0FloatView.put(w0ArrayView[index]);
                 }
             }
         }
@@ -755,10 +991,14 @@ 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++){
-                    v0ArrayView[IX(i,j,k)] = densityArrayView[IX(i,j,k)] * GRAVITY;
+                    index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
+                    u0ArrayView[index] = 0;
+                    v0ArrayView[index] = densityArrayView[index] * GRAVITY;
+                    w0ArrayView[index] = 0;
                 }
             }
         }
@@ -823,6 +1063,10 @@ 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)];
     }
@@ -849,6 +1093,7 @@ 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 = 2;
+        int dim = 10;
         int i = 0;
         long time = 0;
         long lastTime = 0;
@@ -55,7 +55,7 @@ public class Main {
                 //
                 //Simulate
                 //
-                FluidSim.simChunks(simArray,i,0.001f);
+                FluidSim.simChunks(simArray,i,0.01f);
                 time = time + (System.currentTimeMillis() - lastTime);
                 //
                 //Remesh
@@ -70,8 +70,8 @@ public class Main {
                 //redraw
                 GLFWContext.redraw(meshArray);
                 i++;
-                if(i == 1000){
-                    System.out.println(time / 1000.0);
+                if(i == 100){
+                    System.out.println(time / 100.0);
                 }
                 if(i > 3){
                     // scan.next();

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

@@ -2,11 +2,13 @@ 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;
@@ -241,7 +243,7 @@ public class Mesh {
     public static void initShaderProgram(){
         String vsSrc = "";
         ClassLoader classloader = Thread.currentThread().getContextClassLoader();
-        try (BufferedReader is = new BufferedReader(new InputStreamReader(classloader.getResourceAsStream("shader.vs")))){
+        try (BufferedReader is = new BufferedReader(Files.newBufferedReader(new File("C:\\Users\\satellite\\Documents\\fluid-sim\\src\\main\\resources\\shader.vs").toPath()))){
             String temp;
             while((temp = is.readLine())!=null){
                 vsSrc = vsSrc + temp + "\n";
@@ -252,7 +254,7 @@ public class Mesh {
         }
 
         String fsSrc = "";
-        try (BufferedReader is = new BufferedReader(new InputStreamReader(classloader.getResourceAsStream("shader.fs")))){
+        try (BufferedReader is = new BufferedReader(Files.newBufferedReader(new File("C:\\Users\\satellite\\Documents\\fluid-sim\\src\\main\\resources\\shader.fs").toPath()))){
             String temp;
             while((temp = is.readLine())!=null){
                 fsSrc = fsSrc + temp + "\n";