float noise(float2 p,__global float2* rt,__global unsigned char* rti)
{
float2 floor_p = floor(p);//float2(floor(p.x),floor(p.y))
int i=(int)floor_p.x, j=(int)floor_p.y;
float2 n00=rt[hash_idx(i,j,rti)];
float2 n10=rt[hash_idx(i+1,j,rti)];
float2 n01=rt[hash_idx(i,j+1,rti)];
float2 n11=rt[hash_idx(i+1,j+1,rti)];
float2 pf = p - floor_p;
float sx= pow(pf.x,3)*(10-pf.x*(15-pf.x*6));
float sy= pow(pf.y,3)*(10-pf.y*(15-pf.y*6));
return bilerp( dot(pf,n00) , dot(pf-(float2)(1.f,0.f),n10),
dot(pf-(float2)(0.f,1.f),n01) , dot(pf-(float2)(1.f,1.f),n11),
sx, sy);
}
//host
float spin_t = 0.5f*0.03f/0.1f*time_;
for(unsigned int i=0; i<RANDOM_SIZE; ++i){
float theta=rnd_table_spin_[i]*spin_t;
float c=std::cos(theta),s=std::sin(theta);
rnd_table_[i][0]= c*rnd_table_orig_[i][0]+s*rnd_table_orig_[i][1];
rnd_table_[i][1]=-s*rnd_table_orig_[i][0]+c*rnd_table_orig_[i][1];
}
float turb=g*d*0.0026f*noise((px-t*wind_velocity)/0.1f,rt,rti);
^^^^^^^
100k particles
GPU : 10ms
CPU single : 183ms
CPU omp : 93ms
__kernel void update(__global float2* in,
__global float2* rx,
__global float * t,
__global float2* rt,
__global unsigned char* rti){
int const id = get_global_id(0);
float2 mid = in[id];
mid = mid+0.5f*curl_vel(in[id],rx,t[0],rt,rti)*0.04f;
in[id] += curl_vel(mid,rx,t[0],rt,rti)*0.04f;
}float2 curl_vel(float2 x,__global float2* rx,float t,__global float2* rt, __global unsigned char* rti)
{
float const delta_x = (1e-4);
float2 v;
v.x = -(potential((float2)(x.x,x.y+delta_x),rx,t,rt,rti)
- potential((float2)(x.x,x.y-delta_x),rx,t,rt,rti))/(2*delta_x);
v.y = (potential((float2)(x.x+delta_x,x.y),rx,t,rt,rti)
- potential((float2)(x.x-delta_x, x.y),rx,t,rt,rti))/(2*delta_x);
return v;
}float potential(float2 px,__global float2* rx,float t,__global float2* rt, __global unsigned char* rti)
{
float2 wind_velocity = (float2)(0.1f,0.5f);
float const inv_env_sqr = 1/sqrt(0.17f);
float numer= inv_env_sqr*cross2(px,wind_velocity);
float denom= inv_env_sqr;
//---------------hard coded rigid X( ---
float2 dx=px-rx[0];//centrer
float psi=cross2(px,rx[1]);
float phi=((sqrt(dx.x*dx.x)+sqrt(dx.y*dx.y))-0.06f);//radius;
float m=1/(sqrt(phi)+1e-18);
numer+=m*psi;
denom+=m;
//--------------------------------------
float base_psi=numer/denom;
float d=phi/0.1f;
float g=smooth_step(1.1-px.x);
float turb=g*d*0.0026f*noise((px-t*wind_velocity)/0.1f,rt,rti);
return base_psi+turb;
}
//vs
uniform mat4 mvp;
in vec3 in_pos;
in vec2 in_tex;
out vec2 out_tex;
void main(void)
{
gl_Position = mvp * vec4(gl_Vertex.xyz,0.f);
out_tex = in_tex;
}//fs
uniform int mat_size;
uniform sampler2D mat_1;
uniform sampler2D mat_2;
in vec2 texcoord;
void main (void)
{
float result = 0.f;
for(int i = 0;i<mat_size;++i)
{
result += texture2D( mat_1,vec2(i/(float)mat_size,texcoord.y) )
* texture2D( mat_2,vec2(texcoord.y,i/(float)mat_size) );
}
gl_FragColor = vec4(result,1.f);
}//データ作成
GLuint mat_lhs;
::glGenTextures( 1, &mat_lhs );
::glBindTexture( GL_TEXTURE_2D, mat_lhs );
::glTexImage2D(GL_TEXTURE_2D,0,GL_LUMINANCE,mat_size,mat_size,0,GL_LUMINANCE,GL_FLOAT,data);
::glBindTexture( GL_TEXTURE_2D, 0 );//実行
::glEnableClientState(GL_VERTEX_ARRAY);
::glEnableClientState(GL_TEXTURE_COORD_ARRAY);
::glBindBuffer(GL_ARRAY_BUFFER, texed_quad);
::glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, index_buffer);
::glDrawElements(GL_TRIANGLES,buffer_size/sizeof(GLuint),GL_UNSIGNED_INT,0);
::glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
::glBindBuffer(GL_ARRAY_BUFFER, 0);
::glDisableClientState(GL_VERTEX_ARRAY);
::glDisableClientState(GL_TEXTURE_COORD_ARRAY);
::glFlush();
::glReadPixels(0,0,mat_size,mat_size,GL_RGBA,GL_FLOAT,result);
//in raycast function
//*r is ray
int const size = res * center.w;
float4 const mini = center - size*0.5f;
for(int i = 1;1;++i)
{
float4 const p = r.orig_ + r.dir_*i*center.w;
float4 const r_vox_pos = ( p - mini )/center.w;
int4 const r_vox_pos_i = convert_int(r_vox_pos);
//exit without hit
if(r_vox_pos_i.x < 0 || r_vox_pos_i.x >= res-1 ||
r_vox_pos_i.y < 0 || r_vox_pos_i.y >= res-1 ||
r_vox_pos_i.z < 0 || r_vox_pos_i.z >= res)
return false;
sphere sp;
sp.pos_ = convert_float(r_vox_pos_i)*center.w - size/2.f + center.w*0.5f;
sp.r_ = center.w*0.5f;
float t;
if( sphere_intersect(r,sp,&t) )
{
// hit actions
return true;
}