diff --git a/.gitignore b/.gitignore index 77c9be0..8d0637a 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,5 @@ /.settings /.classpath /.project -/.vscode \ No newline at end of file +/.vscode +/shared-folder \ No newline at end of file diff --git a/shared-folder/libfluidsim.dll b/shared-folder/libfluidsim.dll index 2fbdd2c..ae966ff 100644 Binary files a/shared-folder/libfluidsim.dll and b/shared-folder/libfluidsim.dll differ diff --git a/src/main/c/fluidsim.c b/src/main/c/fluidsim.c index 2fdfdba..2ec4474 100644 --- a/src/main/c/fluidsim.c +++ b/src/main/c/fluidsim.c @@ -6,11 +6,11 @@ #define SWAP(x0,x) {float *tmp=x0;x0=x;x=tmp;} #define IX(i,j,k) ((i)+(N)*(j)+(N*N)*(k)) -#define LINEARSOLVERTIMES 20 - +#define LINEARSOLVERTIMES 10 +//https://docs.oracle.com/javase/1.5.0/docs/guide/jni/spec/functions.html void diffuse(int N, int b, float * x, float * x0, float diff, float dt); void advect(int N, int b, float * d, float * d0, float * u, float * v, float * w, float dt); @@ -26,47 +26,43 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate( jint DIM_X, jint DIM_Y, jint DIM_Z, - jfloatArray jx, - jfloatArray jx0, - jfloatArray ju, - jfloatArray jv, - jfloatArray jw, - jfloatArray ju0, - jfloatArray jv0, - jfloatArray jw0, + 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)->GetFloatArrayElements(env,jx,&isCopy); - float * x0 = (*env)->GetFloatArrayElements(env,jx0,&isCopy); - float * u = (*env)->GetFloatArrayElements(env,ju,&isCopy); - float * v = (*env)->GetFloatArrayElements(env,jv,&isCopy); - float * w = (*env)->GetFloatArrayElements(env,jw,&isCopy); - float * u0 = (*env)->GetFloatArrayElements(env,ju0,&isCopy); - float * v0 = (*env)->GetFloatArrayElements(env,jv0,&isCopy); - float * w0 = (*env)->GetFloatArrayElements(env,jw0,&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; + // for ( i=1 ; iSetFloatArrayRegion(env,jx,0,DIM_X*DIM_X*DIM_X,x); - (*env)->SetFloatArrayRegion(env,jx0,0,DIM_X*DIM_X*DIM_X,x0); - (*env)->SetFloatArrayRegion(env,ju,0,DIM_X*DIM_X*DIM_X,u); - (*env)->SetFloatArrayRegion(env,jv,0,DIM_X*DIM_X*DIM_X,v); - (*env)->SetFloatArrayRegion(env,jw,0,DIM_X*DIM_X*DIM_X,w); - (*env)->SetFloatArrayRegion(env,ju0,0,DIM_X*DIM_X*DIM_X,u0); - (*env)->SetFloatArrayRegion(env,jv0,0,DIM_X*DIM_X*DIM_X,v0); - (*env)->SetFloatArrayRegion(env,jw0,0,DIM_X*DIM_X*DIM_X,w0); - // for(int i=0; iSetFloatArrayRegion(env,jx,0,DIM_X*DIM_X*DIM_X,x); - // // jx[IX(i,j,k)] = x[IX(i,j,k)]; - // }}} } - - -int main() { - printf("Hello World from Project Panama Baeldung Article"); - return 0; +void add_source(int N, float * x, float * s, float dt){ + int i; + int size=N*N*N; + for(i=0; iN+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; + for ( i=1 ; iN+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(j0 >= N){ - j0 = N - 1; - } - if(k0 >= N){ - k0 = N - 1; - } - if(i1 >= N){ - i1 = N - 1; - } - if(j1 >= N){ - j1 = N - 1; - } - if(k1 >= N){ - k1 = N - 1; - } - 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)]); - }}} + 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(N, b, d); } + void dens_step(int N, float * x, float * x0, float * u, float * v, float * w, float diff, float dt){ // add_source ( N, x, x0, dt ); - SWAP ( x0, x ); diffuse ( N, 0, x, x0, diff, dt ); - SWAP ( x0, x ); advect ( N, 0, x, x0, u, v, w, dt ); + SWAP(x0, x); + diffuse(N, 0, x, x0, diff, dt); + SWAP(x0, x); + advect(N, 0, x, x0, u, v, w, dt); } void vel_step(int N, float * u, float * v, float * w, float * u0, float * v0, float * w0, float visc, float dt){ - // add_source ( N, u, u0, dt ); add_source ( N, v, v0, dt ); - // SWAP ( u0, u ); diffuse ( N, 1, u, u0, visc, dt ); - // SWAP ( v0, v ); diffuse ( N, 2, v, v0, visc, dt ); - // project ( N, u, v, w, u0, v0 ); - // SWAP ( u0, u ); SWAP ( v0, v ); - // advect ( N, 1, u, u0, u0, v0, dt ); advect ( N, 2, v, v0, u0, v0, dt ); - // project ( N, u, v, w, u0, v0 ); - - // add_source(N, v, x, dt); + add_source(N, u, u0, dt); + add_source(N, v, v0, dt); + add_source(N, w, w0, dt); SWAP(u0, u); diffuse(N, 1, u, u0, visc, dt); SWAP(v0, v); @@ -168,6 +179,7 @@ void project(int N, float * u, float * v, float * w, float * p, float * div){ set_bnd(N, 3, w); } + void lin_solve(int N, int b, float* x, float* x0, float a, float c){ int i, j, k, l; @@ -181,19 +193,6 @@ void lin_solve(int N, int b, float* x, float* x0, float a, float c){ } } -// void set_bnd ( int N, int b, float * x ){ -// int i; -// // for ( i=1 ; i<=N ; i++ ) { -// // x[IX(0 ,i)] = b==1 ? -x[IX(1,i)] : x[IX(1,i)]; -// // x[IX(N+1,i)] = b==1 ? -x[IX(N,i)] : x[IX(N,i)]; -// // x[IX(i,0 )] = b==2 ? -x[IX(i,1)] : x[IX(i,1)]; -// // x[IX(i,N+1)] = b==2 ? -x[IX(i,N)] : x[IX(i,N)]; -// // } -// // x[IX(0 ,0 )] = 0.5*(x[IX(1,0 )]+x[IX(0 ,1)]); -// // x[IX(0 ,N+1)] = 0.5*(x[IX(1,N+1)]+x[IX(0 ,N )]); -// // x[IX(N+1,0 )] = 0.5*(x[IX(N,0 )]+x[IX(N+1,1)]); -// // x[IX(N+1,N+1)] = 0.5*(x[IX(N,N+1)]+x[IX(N+1,N )]); -// } void set_bnd(int N, int b, float * target){ int DIM = N; diff --git a/src/main/c/includes/electrosphere_FluidSim.h b/src/main/c/includes/electrosphere_FluidSim.h index 146de66..5d87a94 100644 --- a/src/main/c/includes/electrosphere_FluidSim.h +++ b/src/main/c/includes/electrosphere_FluidSim.h @@ -10,20 +10,20 @@ extern "C" { #undef electrosphere_FluidSim_DIM #define electrosphere_FluidSim_DIM 18L #undef electrosphere_FluidSim_DIFFUSION_CONSTANT -#define electrosphere_FluidSim_DIFFUSION_CONSTANT 1.0E-4f +#define electrosphere_FluidSim_DIFFUSION_CONSTANT 0.0f #undef electrosphere_FluidSim_VISCOSITY_CONSTANT -#define electrosphere_FluidSim_VISCOSITY_CONSTANT 1.0E-4f +#define electrosphere_FluidSim_VISCOSITY_CONSTANT 0.0f #undef electrosphere_FluidSim_LINEARSOLVERTIMES #define electrosphere_FluidSim_LINEARSOLVERTIMES 20L #undef electrosphere_FluidSim_GRAVITY -#define electrosphere_FluidSim_GRAVITY -10000.0f +#define electrosphere_FluidSim_GRAVITY -1000.0f /* * Class: electrosphere_FluidSim * Method: simulate - * Signature: (III[F[F[F[F[F[F[F[FFFF)V + * Signature: (IIILjava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V */ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate - (JNIEnv *, jobject, jint, jint, jint, jfloatArray, jfloatArray, jfloatArray, jfloatArray, jfloatArray, jfloatArray, jfloatArray, jfloatArray, jfloat, jfloat, jfloat); + (JNIEnv *, jobject, jint, jint, jint, jobject, jobject, jobject, jobject, jobject, jobject, jobject, jobject, jfloat, jfloat, jfloat); #ifdef __cplusplus } diff --git a/src/main/java/electrosphere/FluidSim.java b/src/main/java/electrosphere/FluidSim.java index 969721b..893634c 100644 --- a/src/main/java/electrosphere/FluidSim.java +++ b/src/main/java/electrosphere/FluidSim.java @@ -3,6 +3,8 @@ package electrosphere; import java.io.File; import java.io.IOException; import java.nio.ByteBuffer; +import java.nio.ByteOrder; +import java.nio.FloatBuffer; import java.nio.file.Files; import java.util.LinkedList; import java.util.List; @@ -26,24 +28,50 @@ public class FluidSim { public static final int DIM = 18; - float[] x = new float[DIM * DIM * DIM]; - float[] x0 = new float[DIM * DIM * DIM]; - float[] u = new float[DIM * DIM * DIM]; - float[] v = new float[DIM * DIM * DIM]; - float[] w = new float[DIM * DIM * DIM]; - float[] u0 = new float[DIM * DIM * DIM]; - float[] v0 = new float[DIM * DIM * DIM]; - float[] w0 = new float[DIM * DIM * DIM]; - float[] jacobiAltArray = new float[DIM * DIM * DIM]; + ByteBuffer x = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4); + ByteBuffer x0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4); + ByteBuffer u = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4); + ByteBuffer v = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4); + ByteBuffer w = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4); + ByteBuffer u0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4); + ByteBuffer v0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4); + ByteBuffer w0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4); - static final float DIFFUSION_CONSTANT = 0.0001f; - static final float VISCOSITY_CONSTANT = 0.0001f; + //The densities for every voxel for the current frame + float[] densityArrayView = new float[DIM * DIM * DIM]; + //Should be set to add water to the current frame + float[] density0ArrayView = new float[DIM * DIM * DIM]; + //The force vector for each voxel for the current frame + float[] uArrayView = new float[DIM * DIM * DIM]; + float[] vArrayView = new float[DIM * DIM * DIM]; + 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]; + float[] w0ArrayView = new float[DIM * DIM * DIM]; + + + static final float DIFFUSION_CONSTANT = 0.0f; + static final float VISCOSITY_CONSTANT = 0.0f; static final int LINEARSOLVERTIMES = 20; - static final float GRAVITY = -10000f; + static final float GRAVITY = -1000f; public void setup(){ + x.order(ByteOrder.LITTLE_ENDIAN); + x0.order(ByteOrder.LITTLE_ENDIAN); + u.order(ByteOrder.LITTLE_ENDIAN); + v.order(ByteOrder.LITTLE_ENDIAN); + w.order(ByteOrder.LITTLE_ENDIAN); + u0.order(ByteOrder.LITTLE_ENDIAN); + v0.order(ByteOrder.LITTLE_ENDIAN); + w0.order(ByteOrder.LITTLE_ENDIAN); + FloatBuffer xf = x.asFloatBuffer(); + FloatBuffer uf = u.asFloatBuffer(); + FloatBuffer vf = v.asFloatBuffer(); + FloatBuffer wf = w.asFloatBuffer(); + Random rand = new Random(1); //make a cube of water in the center for(int i = 0; i < DIM; i++){ @@ -54,10 +82,15 @@ public class FluidSim { Math.abs(8 - j) < 4 && Math.abs(10 - k) < 4 ){ - x[IX(i,j,k)] = 1; - u[IX(i,j,k)] = rand.nextFloat() * 0.1f; - v[IX(i,j,k)] = -1f; - w[IX(i,j,k)] = rand.nextFloat() * 0.1f; + xf.put(1); + uf.put(rand.nextFloat() * 0.1f); + vf.put(-1f); + wf.put(rand.nextFloat() * 0.1f); + } else { + xf.put(0); + uf.put(0); + vf.put(0); + wf.put(0); } } } @@ -68,283 +101,167 @@ public class FluidSim { * Runs a frame of the fluid simulation */ public void simulate(int step, float timestep){ - //diffuse vectors - // moveData(oldData[1], data[1]); - // diffuse(timestep, data[1], oldData[1], 1); - // moveData(oldData[2], data[2]); - // diffuse(timestep, data[2], oldData[2], 2); - // moveData(oldData[3], data[3]); - // diffuse(timestep, data[3], oldData[3], 3); - // //project - // project(data[1],data[2],data[3],oldData[1],oldData[2]); - // moveData(oldData[1], data[1]); - // moveData(oldData[2], data[2]); - // moveData(oldData[3], data[3]); - // //advect vectors - // advect(timestep, 1, 1); - // advect(timestep, 2, 2); - // advect(timestep, 3, 3); - // project(data[1],data[2],data[3],oldData[1],oldData[2]); + // + // Add forces and density here + // + addGravity(); - // //diffuse density - // moveData(oldData[0], data[0]); - // diffuse(timestep, data[0], oldData[0], 0); - // //advect density - // moveData(oldData[0], data[0]); - // advect(timestep, 0, 0); + // + //Performs main fluid simulation logic + // + writeNewStateIntoBuffers(); - if(step % 1500 == 0){ - for(int i = 0; i < 4; i++){ - for(int j = 0; j < 4; j++){ - x[IX(15 - i,15,15 - j)] = 1; - v[IX(15 - i,15,15 - j)] = -1; - } - } + if(x.position() > 0){ + x.position(0); + } + if(x0.position() > 0){ + x0.position(0); + } + if(u.position() > 0){ + u.position(0); + } + if(v.position() > 0){ + v.position(0); + } + if(w.position() > 0){ + w.position(0); + } + if(u0.position() > 0){ + u0.position(0); + } + if(v0.position() > 0){ + v0.position(0); + } + if(w0.position() > 0){ + w0.position(0); } - simulate(DIM, DIM, DIM, x, x0, u, v, w, u0, v0, w0, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep); - // vel_step(DIM, DIM, DIM, u, v, w, u0, v0, w0, VISCOSITY_CONSTANT, timestep); + // + //Reads out the results of the fluid sim + // + readDataIntoArrays(); - // dens_step(DIM, DIM, DIM, x, x0, u, v, w, DIFFUSION_CONSTANT, timestep); - - sum(); } - - void lin_solve(int M, int N, int O, int b, float[] x, float[] x0, float a, float c){ - int i, j, k, l; - - // iterate the solver - for ( l=0 ; l 0){ + x.position(0); } - } - - - void diffuse(int M, int N, int O, int b, float[] x, float[] x0, float diff, float dt){ - int max = Math.max(Math.max(M, N), Math.max(N, O)); - float a=dt*diff*max*max*max; - lin_solve(M, N, O, b, x, x0, a, 1+6*a); - } - - void setBounds(float[] target, int b){ - 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)]; - } + if(u.position() > 0){ + u.position(0); } - 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)])); - + if(v.position() > 0){ + v.position(0); } - target[IX(0,0,0)] = (float)((target[IX(1,0,0)]+target[IX(0,1,0)]+target[IX(0,0,1)])/3.0); - target[IX(DIM-1,0,0)] = (float)((target[IX(DIM-2,0,0)]+target[IX(DIM-1,1,0)]+target[IX(DIM-1,0,1)])/3.0); - target[IX(0,DIM-1,0)] = (float)((target[IX(1,DIM-1,0)]+target[IX(0,DIM-2,0)]+target[IX(0,DIM-1,1)])/3.0); - target[IX(0,0,DIM-1)] = (float)((target[IX(0,0,DIM-2)]+target[IX(1,0,DIM-1)]+target[IX(0,1,DIM-1)])/3.0); - target[IX(DIM-1,DIM-1,0)] = (float)((target[IX(DIM-2,DIM-1,0)]+target[IX(DIM-1,DIM-2,0)]+target[IX(DIM-1,DIM-1,1)])/3.0); - target[IX(0,DIM-1,DIM-1)] = (float)((target[IX(1,DIM-1,DIM-1)]+target[IX(0,DIM-2,DIM-1)]+target[IX(0,DIM-1,DIM-2)])/3.0); - target[IX(DIM-1,0,DIM-1)] = (float)((target[IX(DIM-1,0,DIM-2)]+target[IX(DIM-2,0,DIM-1)]+target[IX(DIM-1,1,DIM-1)])/3.0); - target[IX(DIM-1,DIM-1,DIM-1)] = (float)((target[IX(DIM-1,DIM-1,DIM-2)]+target[IX(DIM-1,DIM-2,DIM-1)]+target[IX(DIM-1,DIM-1,DIM-2)])/3.0); - } - - void advect(int M, int N, int O, int b, float[] d, 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*Math.max(Math.max(M, N), Math.max(N, O)); - - for ( i=1 ; iM+0.5f) x=M+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>O+0.5f) z=O+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 >= DIM){ - i0 = DIM - 1; - } - if(j0 >= DIM){ - j0 = DIM - 1; - } - if(k0 >= DIM){ - k0 = DIM - 1; - } - if(i1 >= DIM){ - i1 = DIM - 1; - } - if(j1 >= DIM){ - j1 = DIM - 1; - } - if(k1 >= DIM){ - k1 = DIM - 1; - } - 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)]); - }}} - - setBounds(d, b); - } - - void moveData(float[] set1, float[] set2){ - float tmp = 0; - for(int i = 0; i < set1.length; i++){ - tmp = set1[i]; - set1[i] = set2[i]; - set2[i] = tmp; + if(w.position() > 0){ + w.position(0); } - } - - void project(int M, int N, int O, float[] u, float[] v, float[] w, float[] p, float[] div){ - int i, j, k; - - for ( i=1 ; i 0.0f){ - x[IX(i,j,k)] += dt * GRAVITY * density[IX(i,j,k)]; - } - }}} - } - - void vel_step(int M, int N, int O, float[] u, float[] v, float[] w, float[] u0, float[] v0, float[] w0, float visc, float dt ){ - // add_source ( M, N, O, u, u0, dt ); add_source ( M, N, O, v, v0, dt );add_source ( M, N, O, w, w0, dt ); - add_source(M, N, O, v, x, dt); - moveData(u0, u); - diffuse(M, N, O, 1, u, u0, visc, dt); - moveData(v0, v); - diffuse(M, N, O, 2, v, v0, visc, dt); - moveData(w0, w); - diffuse(M, N, O, 3, w, w0, visc, dt); - project(M, N, O, u, v, w, u0, v0); - moveData(u0, u); - moveData(v0, v); - moveData(w0, w); - advect(M, N, O, 1, u, u0, u0, v0, w0, dt); - advect(M, N, O, 2, v, v0, u0, v0, w0, dt); - advect(M, N, O, 3, w, w0, u0, v0, w0, dt); - project(M, N, O, u, v, w, u0, v0); - } - - void sum(){ - double[] sum = new double[8]; for(int i = 0; i < DIM; i++){ for(int j = 0; j < DIM; j++){ for(int k = 0; k < DIM; k++){ - sum[0] = sum[0] + x[IX(i,j,k)]; - sum[1] = sum[1] + u[IX(i,j,k)]; - sum[2] = sum[2] + v[IX(i,j,k)]; - sum[3] = sum[3] + w[IX(i,j,k)]; + densityArrayView[IX(i,j,k)] = x.getFloat(); + uArrayView[IX(i,j,k)] = u.getFloat(); + vArrayView[IX(i,j,k)] = v.getFloat(); + wArrayView[IX(i,j,k)] = w.getFloat(); } } } - // System.out.println(sum[0] + " " + sum[1] + " " + sum[2] + " " + sum[3]); } + /** + * Writes data from the java-side arrays into buffers that get passed into c-side + */ + private void writeNewStateIntoBuffers(){ + if(x0.position() > 0){ + x0.position(0); + } + if(u0.position() > 0){ + u0.position(0); + } + if(v0.position() > 0){ + v0.position(0); + } + if(w0.position() > 0){ + w0.position(0); + } + FloatBuffer x0FloatView = x0.asFloatBuffer(); + FloatBuffer u0FloatView = u0.asFloatBuffer(); + FloatBuffer v0FloatView = v0.asFloatBuffer(); + FloatBuffer w0FloatView = w0.asFloatBuffer(); + 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)]); + } + } + } + } + + /** + * Adds gravity to the simulation + */ + private void addGravity(){ + 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; + } + } + } + } + + /** + * The native function call to simulate a frame of fluid + * @param DIM_X + * @param DIM_Y + * @param DIM_Z + * @param x + * @param x0 + * @param u + * @param v + * @param w + * @param u0 + * @param v0 + * @param w0 + * @param DIFFUSION_CONSTANT + * @param VISCOSITY_CONSTANT + * @param timestep + */ private native void simulate( int DIM_X, int DIM_Y, int DIM_Z, - float[] x, - float[] x0, - float[] u, - float[] v, - float[] w, - float[] u0, - float[] v0, - float[] w0, + ByteBuffer x, + ByteBuffer x0, + ByteBuffer u, + ByteBuffer v, + ByteBuffer w, + ByteBuffer u0, + ByteBuffer v0, + ByteBuffer w0, float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep ); public float[] getData(){ - return x; + return densityArrayView; } public static final int IX(int x, int y, int z){ return ((x)+(DIM)*(y) + (DIM)*(DIM)*(z)); } - static final int[] vectorOffsetMap = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)]; - static final int[] vectorOffsetMapXN = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)]; - static final int[] vectorOffsetMapXP = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)]; - static final int[] vectorOffsetMapYN = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)]; - static final int[] vectorOffsetMapYP = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)]; - static final int[] vectorOffsetMapZN = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)]; - static final int[] vectorOffsetMapZP = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)]; - static { - int i = 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++){ - // if(x != y || y != z){ - vectorOffsetMap[i] = IX(x,y,z); - vectorOffsetMapXN[i] = IX(x-1, y, z); - vectorOffsetMapXP[i] = IX(x+1, y, z); - vectorOffsetMapYN[i] = IX( x,y-1, z); - vectorOffsetMapYP[i] = IX( x,y+1, z); - vectorOffsetMapZN[i] = IX( x, y,z-1); - vectorOffsetMapZP[i] = IX( x, y,z+1); - i++; - // } - } - } - } - } - }