ProjectPhysX
commented
11 months ago Hi @ma-sadeghi,
these are super old setups from one of the first prototype versions. I managed to find the old setup source code in one of my backups. Be aware it's very messy code and not out-of-the-box compatible with the clean new version. Nontheless I hope you can salvage what you are looking for!
Kind regards, Moritz
void initializeFluid() {
/*// random positions
const int N = 300;
float** p = new float*[N];
for(int i=0; i<N; i++) {
p[i] = new float[3];
p[i][0] = random(lattice.size_x);
p[i][1] = random(lattice.size_y);
p[i][2] = random(lattice.size_z);
}/**/
const uint sx=lattice.size_x, sy=lattice.size_y, sz=lattice.size_z, s=lattice.point_number;
for(uint x=0; x<sx; x++) {
for(uint y=0; y<sy; y++) {
for(uint z=0; z<sz; z++) {
const uint n = x+(y+z*sy)*sx; // set u[ n], u[ s+n], u[2*s+n], rho[n], flags[n]
#ifndef VALIDATION
#ifndef D2Q9
// box
//if(x==0||x==sx-1||y==0||y==sy-1||z==0||z==sz-1) flags[n] = 1; // all non periodic
//if(x==0||x==sx-1) flags[n] = 1; // x non periodic
//if(y==0||y==sy-1) flags[n] = 1; // y non periodic
//if(z==0||z==sz-1) flags[n] = 1; // z non periodic
// pressure gradient
/*if(y== 0) rho[n] = 1.3f;
if(y==sy-1) rho[n] = 0.7f;
if(y== 0) flags[n] = 2;
if(y==sy-1) flags[n] = 2;*/
// lid-driven cavity
if(z==sz-1) u[ s+n] = 0.5f;
if(y== 0) flags[n] = 1; // 3D
if(y==sy-1) flags[n] = 1;
if(z== 0) flags[n] = 1;
if(x==sx-1) flags[n] = 1;
if(x== 0) flags[n] = 1;
if(z==sz-1) flags[n] = 2;/**/
/*if(y==sy/2&&(z<sz/8&&z>sz/16)) flags[n] = 1;
/if(y== 0) u[n] = 0.5f;
if(x== 0) flags[n] = 1; // 2D
if(x==sx-1) flags[n] = 1;
if(y==sy-1) flags[n] = 1;
if(y== 0) flags[n] = 2;*/
// Taylor couette flow
/*if(!barZ((float)x, (float)y, (float)z, (float)sx/2, (float)sy/2, (float)sz/2, (float)sx/2-1, (float)sz)) flags[n] = 1;
if( barZ((float)x, (float)y, (float)z, (float)sx/2, (float)sy/2, (float)sz/2, (float)sx/4, (float)sz)) flags[n] = 1;
if(pipeZ((float)x, (float)y, (float)z, (float)sx/2, (float)sy/2, (float)sz/2, (float)sx/4, (float)sz)) {
u[ n] = ((float)y-(float)sy/2)*0.01f;
u[ s+n] = -((float)x-(float)sx/2)*0.01f;
u[2*s+n] = (1.0f-random(2.0f))*0.001f;
flags[n] = 2;
}/**/
// pipes
/*const int N = 3;
const float D = (float)max(sx, sz);
const float L = (float)sz;
const float r = D/(2.0f*N);
const float l = sy-3.0f*2.0f*r;
const float py = 0.5f*L;
for(int j=0; j<N; j++) {
const float pz = j*D/N+r;
for(int i=0; i<N; i++) {
const float px = i*D/N+r;
if(pipeY(x, y, z, px, py, pz, r, r, L-4*r+1)) flags[n] = 1;
if(spiralY(x, y, z, px, py, pz, r, r, L-4*r+1)) flags[n] = 1;
if(i!=N-1 && i%2==0 && y<=2*r && torusZ(x, y, z, px+r, 2*r, pz, r, r)) flags[n] = 1;
if(i!=0 && i%2==0 && y>=L-2*r && torusZ(x, y, z, px-r, L-2*r, pz, r, r)) flags[n] = 1;
}
if(j!=N-1 && j%2==0 && y<=2*r && torusX(x, y, z, L-r, 2*r, pz+r, r, r)) flags[n] = 1;
if(j!=0 && j%2==0 && y>=L-2*r && torusX(x, y, z, r, L-2*r, pz-r, r, r)) flags[n] = 1;
}
if(y== 0) {
if(sq(x-(L-r))+sq(z-(L-r))<sq(r)) {
rho[n] = 3.0f;
flags[n] = 2;
} else {
flags[n] = 1;
}
}
if(y==sy-1) {
if(sq(x-r)+sq(z-r)<sq(r)) {
rho[n] = 0.9f;
flags[n] = 2;
} else {
flags[n] = 1;
}
}
if(pipeY(x, y, z, r, L-r, r, r, r, 2*r)) flags[n] = 1;
if(pipeY(x, y, z, L-r, r, L-r, r, r, 2*r)) flags[n] = 1;
if(spiralY(x, y, z, r, L-r, r, r, r, 2*r)) flags[n] = 1;
if(spiralY(x, y, z, L-r, r, L-r, r, r, 2*r)) flags[n] = 1;/**/
// porous medium
/*if(x==0||x==sx-1) flags[n] = 1; // x non periodic
if(z==0||z==sz-1) flags[n] = 1; // z non periodic
for(int i=0; i<N; i++) {
if(sphere(x, y, z, p[i][0], p[i][1], p[i][2], 10.0f)) {
flags[n] = 1;
break;
}
}/**/
//if(sq(x-sx/2)+sq(z-sz/2)>=sq(min(sz, sz)/2-1)) flags[n] = 1; // round channel
//if(sphere(x, y, z, sx/2, sy/2, sz/2, 38.0f)) flags[n] = 1; // sphere turbulence
// temperature
/*if(z==1 && sq(x-sx/4)+sq(y-sy/4)<sq(20)) {
T[n] = 1.9f;
flags[n] = 3;
}
if(z==sz-2) {
T[n] = 0.9f;
flags[n] = 3;
}
if(x==0||x==sx-1||y==0||y==sy-1||z==0||z==sz-1) flags[n] = 1;*/
// heat
/*if(z==1) {
flags[n] = 2;
float r = sqrt(sq(x-sx/2)+sq(y-sy/2)+0.01);
u[ n] = -0.32f*(x-sx/2)/r;
u[ s+n] = -0.32f*(y-sy/2)/r;
}*/
// inflow/outflow
/*if(x==0||x==sx-1||y==0||y==sy-1||z==0||z==sz-1) flags[n] = 1; // all non periodic
if(x==0&&y>sy*0.75f&&y<sy*0.8f&&z>sz*0.1f&&z<sz*0.15f) u[ n] = 0.6f;
if(y==0&&x>sx*0.2f&&x<sx*0.25f&&z>sz*0.1f&&z<sz*0.15f) u[ s+n] = 0.6f;
if(z==sz-1&&x>sx*0.75f&&x<sx*0.8f&&y>sy*0.75f&&y<sy*0.8f) u[2*s+n] = -0.6f;
if(y==sy-1 && sq(x-sx*0.5f)+sq(z-sz*0.5f)<sq(3.0f)) u[ s+n] = -0.6f;
//if(y==sy-1 && sq(x-sx*0.5f)+sq(z-sz*0.5f)>sq(4.0f) && sq(x-sx*0.5f)+sq(z-sz*0.5f)<sq(7.0f)) u[ s+n] = -0.4f;
if(x==0&&y>sy*0.75f&&y<sy*0.8f&&z>sz*0.1f&&z<sz*0.15f) flags[n] = 2;
if(y==0&&x>sx*0.2f&&x<sx*0.25f&&z>sz*0.1f&&z<sz*0.15f) flags[n] = 2;
if(z==sz-1&&x>sx*0.75f&&x<sx*0.8f&&y>sy*0.75f&&y<sy*0.8f) flags[n] = 2;
if(y==sy-1 && sq(x-sx*0.5f)+sq(z-sz*0.5f)<sq(3.0f)) flags[n] = 2;
//if(y==sy-1 && sq(x-sx*0.5f)+sq(z-sz*0.5f)>sq(4.0f) && sq(x-sx*0.5f)+sq(z-sz*0.5f)<sq(7.0f)) flags[n] = 2;*/
// rocket nozzle
/*u[2*s+n] = -0.1f;
rho[n] = 0.9f;
if(z==sz-1) rho[n] = 1.0f;
if(z== 0) rho[n] = 0.6f;
if(z==sz-2 && sq(x-sx*0.5f)+sq(y-sy*0.5f)<sq(6)) {
u[2*s+n] = -0.1f;
rho[n] = 30.0f;
}
//if(z==sz-1) flags[n] = 2;
if(z==sz-2 && sq(x-sx*0.5f)+sq(y-sy*0.5f)<sq(7)) flags[n] = 2;
if(z==sz-1 && sq(x-sx*0.5f)+sq(y-sy*0.5f)<sq(8)) flags[n] = 1;
if(z<=sz-1 && z>=sz-1-50 && sq(x-sx*0.5f)+sq(y-sy*0.5f)<=sq(3.2f*sqrtf(sz-1-z)+4+2.3f) && sq(x-sx*0.5f)+sq(y-sy*0.5f)>=sq(3.2f*sqrtf(sz-1-z)+4)) flags[n] = 1;
if(z==sz-1) flags[n] = 1;
if(z==0) flags[n] = 2;/**/
// three pipes
/*if(pipeX(x, y, z, sx*0.5f, 30, sz*0.5f, 15, 48)) flags[n] = 1;
if(pipeX(x, y, z, sx*0.5f, 80, sz*0.5f-20, 15, 48)) flags[n] = 1;
if(pipeX(x, y, z, sx*0.5f, 65, sz*0.5f+20, 15, 48)) flags[n] = 1;*/
// relax 1
/*u[ n] = 0.001f*(sz-z)*(y-sy*0.5f);
u[ s+n] = -0.001f*(sz-z)*(x-sx*0.5f);
u[2*s+n] = 0.0f;*/
// relax 2
/*float v = 1.0f;
float r = sqrtf(sq(x-sx*0.5f)+sq(y-sy*0.5f)+sq(z-sz*0.5f))+1.0f;
u[ n] = v*(z-sz)/sz*(y-sy*0.5f)/r;
u[ s+n] = -v*(z-sz)/sz*(x-sx*0.5f)/r;
u[2*s+n] = 0.0f;*/
// relax 3
/*u[ n] = 0.001f*(sz-z)*(y-0.3f*sy);
u[ s+n] = -0.001f*(sz-z)*(x-0.3f*sx);
u[2*s+n] = 0.0f;*/
// cones
/*if(sphere(x, y, z, sx*0.5f, sy*0.5f, sz*0.5f, 40) && y<sy*0.5f) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, 20, sz*0.5f, 70, 48, 40)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, 60, sz*0.5f, 48, 40)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, 100, sz*0.5f, 48, 20, 40)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f, sz*0.5f, 20, sy-240)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy-100, sz*0.5f, 20, 48, 40)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy-60, sz*0.5f, 48, 40)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy-20, sz*0.5f, 48, 70, 40)) flags[n] = 1;
//if(pipeY(x, y, z, sx*0.5f, sy*0.5f-60, sz*0.5f, 0, 20, 50)) flags[n] = 1;
//if(pipeY(x, y, z, sx*0.5f, sy*0.5f, sz*0.5f, 20, 70)) flags[n] = 1;
//if(pipeY(x, y, z, sx*0.5f, sy*0.5f+60, sz*0.5f, 20, 0, 50)) flags[n] = 1;
//if(pipeY(x, y, z, sx*0.5f, sy*0.5f-115, sz*0.5f, 80, 30, 100)) flags[n] = 1;
//if(pipeY(x, y, z, sx*0.5f, sy*0.5f-15, sz*0.5f, 30, 80, 100)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f-20, sz*0.5f, 25, 40, 50)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f- 0, sz*0.5f, 20, 35, 50)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f+20, sz*0.5f, 15, 30, 50)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f+40, sz*0.5f, 10, 25, 50)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f+60, sz*0.5f, 5, 20, 50)) flags[n] = 1;
//if(pipeY(x, y, z, sx*0.5f, sy*0.5f+115, sz*0.5f, 30, 80, 100)) flags[n] = 1;
//if(pipeY(x, y, z, sx*0.5f, sy*0.5f+15, sz*0.5f, 80, 30, 100)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f+20, sz*0.5f, 40, 25, 50)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f+ 0, sz*0.5f, 35, 20, 50)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f-20, sz*0.5f, 30, 15, 50)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f-40, sz*0.5f, 25, 10, 50)) flags[n] = 1;
if(pipeY(x, y, z, sx*0.5f, sy*0.5f-60, sz*0.5f, 20, 5, 50)) flags[n] = 1;
//if(y>0&&y%20==0 && sq(x-sx*0.5f)+sq(z-sz*0.5f)<sq(40)) flags[n] = 1;
//if(y%20==10 && sq(x-sx*0.5f)+sq(z-sz*0.5f)>sq(25)) flags[n] = 1;*/
// free surface
#ifdef SURFACE
/*if(x<sx/2&&y<sy/2) {
flags[n] = TYPE_F;
//u[2*s+n] = -0.1f;
} else if((x<sx/2&&y==sy/2)||(x==sx/2&&y<sy/2)) {
flags[n] = TYPE_I;
//u[2*s+n] = -0.1f;
} else {
flags[n] = TYPE_G;
}*/
const float radius = 48;
const float sqr = sq(x-sx*0.5f)+sq(y-sy*0.5f-40)+sq(z-sz*0.5f);
if(sqr<sq(radius)) {
flags[n] |= TYPE_F;
u[ s+n] = -0.3f;
u[2*s+n] = 0.3f;
} else if(sqr>sq(radius)&&sqr<sq(radius+1)) {
flags[n] |= TYPE_I;
u[ s+n] = -0.3f;
u[2*s+n] = 0.3f;
} else {
flags[n] |= TYPE_G;
}
if(x==0||x==sx-1) flags[n] = TYPE_W; // x non periodic
if(y==0||y==sy-1) flags[n] = TYPE_W; // y non periodic
if(z==0||z==sz-1) flags[n] = TYPE_W; // z non periodic
#endif
#else
// 2D
u[n] = 0.1f;
if(x== 0) rho[n] = 1.2f;
if(x==sx-1) rho[n] = 0.8f;
if(x== 0) flags[n] = 2;
if(x==sx-1) flags[n] = 2;
if(y== 0) flags[n] = 1;
if(y==sy-1) flags[n] = 1;
#endif
#else
#ifndef D2Q9
//u[s+n] = 0.25f;
if(sq(x-sx*0.5f+0.5f)+sq(z-sz*0.5f+0.5f)>=sq(min(sx,sz)*0.5f-1.0f)) flags[n] = 1; // 3D
#else
//u[n] = 0.25f;
if(y==0||y==sy-1) flags[n] = 1; // 2D
#endif
#endif
}
}
}
/*// random positions
for(int i=0; i<N; i++) {
delete[] p[i];
}
delete[] p;/**/
}
Thanks for creating this amazing project. I was wondering if you'd be willing to share some of the example scripts that you've showcased on the YouTube channel. Specifically, the porous medium flow in this video (demo 5). Thanks!