studiorailgun/fluid-sim
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

C code works!

Browse Source
This commit is contained in:
unknown 2023-07-16 10:14:31 -04:00
parent 63e3eeba79
commit 8e6b6e4095
5 changed files with 270 additions and 353 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@ -5,4 +5,5 @@
/.settings /.settings
/.classpath /.classpath
/.project /.project
/.vscode /.vscode
/shared-folder

BIN
shared-folder/libfluidsim.dll
View File

Binary file not shown.

177
src/main/c/fluidsim.c
View File

@ -6,11 +6,11 @@
#define SWAP(x0,x) {float *tmp=x0;x0=x;x=tmp;} #define SWAP(x0,x) {float *tmp=x0;x0=x;x=tmp;}
#define IX(i,j,k) ((i)+(N)*(j)+(N*N)*(k)) #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 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); 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_X,
jint DIM_Y, jint DIM_Y,
jint DIM_Z, jint DIM_Z,
jfloatArray jx, jobject jx,
jfloatArray jx0, jobject jx0,
jfloatArray ju, jobject ju,
jfloatArray jv, jobject jv,
jfloatArray jw, jobject jw,
jfloatArray ju0, jobject ju0,
jfloatArray jv0, jobject jv0,
jfloatArray jw0, jobject jw0,
jfloat DIFFUSION_RATE, jfloat DIFFUSION_RATE,
jfloat VISCOSITY_RATE, jfloat VISCOSITY_RATE,
jfloat timestep){ jfloat timestep){
jboolean isCopy; jboolean isCopy;
float * x = (*env)->GetFloatArrayElements(env,jx,&isCopy); float * x = (*env)->GetDirectBufferAddress(env,jx);
float * x0 = (*env)->GetFloatArrayElements(env,jx0,&isCopy); float * x0 = (*env)->GetDirectBufferAddress(env,jx0);
float * u = (*env)->GetFloatArrayElements(env,ju,&isCopy); float * u = (*env)->GetDirectBufferAddress(env,ju);
float * v = (*env)->GetFloatArrayElements(env,jv,&isCopy); float * v = (*env)->GetDirectBufferAddress(env,jv);
float * w = (*env)->GetFloatArrayElements(env,jw,&isCopy); float * w = (*env)->GetDirectBufferAddress(env,jw);
float * u0 = (*env)->GetFloatArrayElements(env,ju0,&isCopy); float * u0 = (*env)->GetDirectBufferAddress(env,ju0);
float * v0 = (*env)->GetFloatArrayElements(env,jv0,&isCopy); float * v0 = (*env)->GetDirectBufferAddress(env,jv0);
float * w0 = (*env)->GetFloatArrayElements(env,jw0,&isCopy); float * w0 = (*env)->GetDirectBufferAddress(env,jw0);
int N = DIM_X;
int i,j,k;
// for ( i=1 ; i<N-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<N-1 ; k++ ) {
// if(v[IX(i,j,k)] < -0.5){
// printf("%d %d %d %f \n",i,j,k,v[IX(i,j,k)]);
// }
// }}}
vel_step(DIM_X, u, v, w, u0, v0, w0, VISCOSITY_RATE, timestep); vel_step(DIM_X, u, v, w, u0, v0, w0, VISCOSITY_RATE, timestep);
dens_step(DIM_X, x, x0, u, v, w, DIFFUSION_RATE, timestep); dens_step(DIM_X, x, x0, u, v, w, DIFFUSION_RATE, timestep);
(*env)->SetFloatArrayRegion(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; i<DIM_X; i++){for(int j=0; j<DIM_X; j++){for(int k=0; k<DIM_X; k++){
// (*env)->SetFloatArrayRegion(env,jx,0,DIM_X*DIM_X*DIM_X,x);
// // jx[IX(i,j,k)] = x[IX(i,j,k)];
// }}}
} }
void add_source(int N, float * x, float * s, float dt){
int i;
int main() { int size=N*N*N;
printf("Hello World from Project Panama Baeldung Article"); for(i=0; i<size; i++){
return 0; x[i] += dt*s[i];
}
} }
@ -77,57 +73,72 @@ 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){ void advect(int N, int b, float * d, float * d0, float * u, float * v, float * w, float dt){
int i, j, k, i0, j0, k0, i1, j1, k1; int i, j, k, i0, j0, k0, i1, j1, k1;
float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz; float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz;
dtx=dty=dtz=dt*N; dtx=dty=dtz=dt*N;
for ( i=1 ; i<N-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<N-1 ; k++ ) { for ( i=1 ; i<N-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<N-1 ; k++ ) {
x = i-dtx*u[IX(i,j,k)]; y = j-dty*v[IX(i,j,k)]; z = k-dtz*w[IX(i,j,k)]; 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 (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 (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; 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; s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1; u1 = z-k0; u0 = 1-u1;
if(i0 >= N){ if(i0 >= N){
i0 = N - 1; i0 = N - 1;
} }
if(j0 >= N){ // if(i0 < 0){
j0 = N - 1; // i0 = 0;
} // }
if(k0 >= N){ if(j0 >= N){
k0 = N - 1; j0 = N - 1;
} }
if(i1 >= N){ // if(j0 < 0){
i1 = N - 1; // j0 = 0;
} // }
if(j1 >= N){ if(k0 >= N){
j1 = N - 1; k0 = N - 1;
} }
if(k1 >= N){ // if(k0 < 0){
k1 = N - 1; // k0 = 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)])+ if(i1 >= N){
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)]); 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); set_bnd(N, b, d);
} }
void dens_step(int N, float * x, float * x0, float * u, float * v, float * w, float diff, float dt){ 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 ); // add_source ( N, x, x0, dt );
SWAP ( x0, x ); diffuse ( N, 0, x, x0, diff, dt ); SWAP(x0, x);
SWAP ( x0, x ); advect ( N, 0, x, x0, u, v, w, dt ); 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){ 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 ); add_source(N, u, u0, dt);
// SWAP ( u0, u ); diffuse ( N, 1, u, u0, visc, dt ); add_source(N, v, v0, dt);
// SWAP ( v0, v ); diffuse ( N, 2, v, v0, visc, dt ); add_source(N, w, w0, 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);
SWAP(u0, u); SWAP(u0, u);
diffuse(N, 1, u, u0, visc, dt); diffuse(N, 1, u, u0, visc, dt);
SWAP(v0, v); 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); set_bnd(N, 3, w);
} }
void lin_solve(int N, int b, float* x, float* x0, float a, float c){ void lin_solve(int N, int b, float* x, float* x0, float a, float c){
int i, j, k, l; 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){ void set_bnd(int N, int b, float * target){
int DIM = N; int DIM = N;

10
src/main/c/includes/electrosphere_FluidSim.h
View File

@ -10,20 +10,20 @@ extern "C" {
#undef electrosphere_FluidSim_DIM #undef electrosphere_FluidSim_DIM
#define electrosphere_FluidSim_DIM 18L #define electrosphere_FluidSim_DIM 18L
#undef electrosphere_FluidSim_DIFFUSION_CONSTANT #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 #undef electrosphere_FluidSim_VISCOSITY_CONSTANT
#define electrosphere_FluidSim_VISCOSITY_CONSTANT 1.0E-4f #define electrosphere_FluidSim_VISCOSITY_CONSTANT 0.0f
#undef electrosphere_FluidSim_LINEARSOLVERTIMES #undef electrosphere_FluidSim_LINEARSOLVERTIMES
#define electrosphere_FluidSim_LINEARSOLVERTIMES 20L #define electrosphere_FluidSim_LINEARSOLVERTIMES 20L
#undef electrosphere_FluidSim_GRAVITY #undef electrosphere_FluidSim_GRAVITY
#define electrosphere_FluidSim_GRAVITY -10000.0f #define electrosphere_FluidSim_GRAVITY -1000.0f
/* /*
* Class: electrosphere_FluidSim * Class: electrosphere_FluidSim
* Method: simulate * 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 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 #ifdef __cplusplus
} }

433
src/main/java/electrosphere/FluidSim.java
View File

@ -3,6 +3,8 @@ package electrosphere;
import java.io.File; import java.io.File;
import java.io.IOException; import java.io.IOException;
import java.nio.ByteBuffer; import java.nio.ByteBuffer;
import java.nio.ByteOrder;
import java.nio.FloatBuffer;
import java.nio.file.Files; import java.nio.file.Files;
import java.util.LinkedList; import java.util.LinkedList;
import java.util.List; import java.util.List;
@ -26,24 +28,50 @@ public class FluidSim {
public static final int DIM = 18; public static final int DIM = 18;
float[] x = new float[DIM * DIM * DIM]; ByteBuffer x = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
float[] x0 = new float[DIM * DIM * DIM]; ByteBuffer x0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
float[] u = new float[DIM * DIM * DIM]; ByteBuffer u = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
float[] v = new float[DIM * DIM * DIM]; ByteBuffer v = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
float[] w = new float[DIM * DIM * DIM]; ByteBuffer w = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
float[] u0 = new float[DIM * DIM * DIM]; ByteBuffer u0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
float[] v0 = new float[DIM * DIM * DIM]; ByteBuffer v0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
float[] w0 = new float[DIM * DIM * DIM]; ByteBuffer w0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
float[] jacobiAltArray = new float[DIM * DIM * DIM];
static final float DIFFUSION_CONSTANT = 0.0001f; //The densities for every voxel for the current frame
static final float VISCOSITY_CONSTANT = 0.0001f; 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 int LINEARSOLVERTIMES = 20;
static final float GRAVITY = -10000f; static final float GRAVITY = -1000f;
public void setup(){ 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); Random rand = new Random(1);
//make a cube of water in the center //make a cube of water in the center
for(int i = 0; i < DIM; i++){ for(int i = 0; i < DIM; i++){
@ -54,10 +82,15 @@ public class FluidSim {
Math.abs(8 - j) < 4 && Math.abs(8 - j) < 4 &&
Math.abs(10 - k) < 4 Math.abs(10 - k) < 4
){ ){
x[IX(i,j,k)] = 1; xf.put(1);
u[IX(i,j,k)] = rand.nextFloat() * 0.1f; uf.put(rand.nextFloat() * 0.1f);
v[IX(i,j,k)] = -1f; vf.put(-1f);
w[IX(i,j,k)] = rand.nextFloat() * 0.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 * Runs a frame of the fluid simulation
*/ */
public void simulate(int step, float timestep){ 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]); //Performs main fluid simulation logic
// diffuse(timestep, data[0], oldData[0], 0); //
// //advect density writeNewStateIntoBuffers();
// moveData(oldData[0], data[0]);
// advect(timestep, 0, 0);
if(step % 1500 == 0){ if(x.position() > 0){
for(int i = 0; i < 4; i++){ x.position(0);
for(int j = 0; j < 4; j++){ }
x[IX(15 - i,15,15 - j)] = 1; if(x0.position() > 0){
v[IX(15 - i,15,15 - j)] = -1; 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); 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; * Used post-simulation call to read data from the simulation from buffers back into arrays
*/
// iterate the solver private void readDataIntoArrays(){
for ( l=0 ; l<LINEARSOLVERTIMES ; l++ ) { //
// update for each cell //Read data back into arrays
for ( i=1 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) { //
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; if(x.position() > 0){
}}} x.position(0);
setBounds(x, b);
} }
} if(u.position() > 0){
u.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)];
}
} }
for(int x = 1; x < DIM-1; x++){ if(v.position() > 0){
target[IX(x,0,0)] = (float)(0.5f * (target[IX(x,1,0)] + target[IX(x,0,1)])); v.position(0);
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); if(w.position() > 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); w.position(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 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
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>M+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;
} }
}
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<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
div[IX(i,j,k)] = (float)(-1.0/3.0*((u[IX(i+1,j,k)]-u[IX(i-1,j,k)])/M+(v[IX(i,j+1,k)]-v[IX(i,j-1,k)])/M+(w[IX(i,j,k+1)]-w[IX(i,j,k-1)])/M));
p[IX(i,j,k)] = 0;
}}}
setBounds(div, 0);
setBounds(p, 0);
lin_solve(M, N, O, 0, p, div, 1, 6);
for ( i=1 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
u[IX(i,j,k)] -= 0.5f*M*(p[IX(i+1,j,k)]-p[IX(i-1,j,k)]);
v[IX(i,j,k)] -= 0.5f*M*(p[IX(i,j+1,k)]-p[IX(i,j-1,k)]);
w[IX(i,j,k)] -= 0.5f*M*(p[IX(i,j,k+1)]-p[IX(i,j,k-1)]);
}}}
setBounds(u, 1);
setBounds(v, 2);
setBounds(w, 3);
}
void dens_step(int M, int N, int O, float[] x, float[] x0, float[] u, float[] v, float[] w, float diff, float dt){
// add_source ( M, N, O, x, x0, dt );
moveData ( x0, x ); diffuse ( M, N, O, 0, x, x0, diff, dt );
moveData ( x0, x ); advect ( M, N, O, 0, x, x0, u, v, w, dt );
}
void add_source(int M, int N, int O, float[] x, float[] density, float dt){
int i, j, k;
for ( i=1 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
if(density[IX(i,j,k)] > 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 i = 0; i < DIM; i++){
for(int j = 0; j < DIM; j++){ for(int j = 0; j < DIM; j++){
for(int k = 0; k < DIM; k++){ for(int k = 0; k < DIM; k++){
sum[0] = sum[0] + x[IX(i,j,k)]; densityArrayView[IX(i,j,k)] = x.getFloat();
sum[1] = sum[1] + u[IX(i,j,k)]; uArrayView[IX(i,j,k)] = u.getFloat();
sum[2] = sum[2] + v[IX(i,j,k)]; vArrayView[IX(i,j,k)] = v.getFloat();
sum[3] = sum[3] + w[IX(i,j,k)]; 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( private native void simulate(
int DIM_X, int DIM_X,
int DIM_Y, int DIM_Y,
int DIM_Z, int DIM_Z,
float[] x, ByteBuffer x,
float[] x0, ByteBuffer x0,
float[] u, ByteBuffer u,
float[] v, ByteBuffer v,
float[] w, ByteBuffer w,
float[] u0, ByteBuffer u0,
float[] v0, ByteBuffer v0,
float[] w0, ByteBuffer w0,
float DIFFUSION_CONSTANT, float DIFFUSION_CONSTANT,
float VISCOSITY_CONSTANT, float VISCOSITY_CONSTANT,
float timestep float timestep
); );
public float[] getData(){ public float[] getData(){
return x; return densityArrayView;
} }
public static final int IX(int x, int y, int z){ public static final int IX(int x, int y, int z){
return ((x)+(DIM)*(y) + (DIM)*(DIM)*(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++;
// }
}
}
}
}
} }