Gonkee's Epic Physics Engine 3D (Gepe3D)
Computational fluid dynamics + rigid bodies + soft bodies, made in C# + OpenGL + OpenCL
Click here to download ZIP file containing EXE
(or click the releases section on the right)
public void Update(float delta, KeyboardState keyboardState) {
Vector3 camOffset = new Vector3(7, 0, 0);
camOffset = Vector3.TransformColumn( Matrix3.CreateRotationZ( MathHelper.DegreesToRadians(pitch) ), camOffset );
camOffset = Vector3.TransformColumn( Matrix3.CreateRotationY( MathHelper.DegreesToRadians(yaw) ), camOffset );
centrePos = particleSystem.GetPos(centreParticleID);
activeCam.SetPos(centrePos + camOffset);
activeCam.LookAt(centrePos);
activeCam.UpdateLocalVectors();
Vector3 movement = new Vector3();
if ( keyboardState.IsKeyDown(Keys.W) ) { movement.X += 1; }
if ( keyboardState.IsKeyDown(Keys.S) ) { movement.X -= 1; }
if ( keyboardState.IsKeyDown(Keys.A) ) { movement.Z -= 1; }
if ( keyboardState.IsKeyDown(Keys.D) ) { movement.Z += 1; }
if ( keyboardState.IsKeyDown(Keys.Space) ) { movement.Y += 1; }
if ( keyboardState.IsKeyDown(Keys.LeftShift) ) { movement.Y -= 1; }
if (movement.Length != 0) movement.Normalize();
Matrix4 rotation = Matrix4.CreateRotationY( MathHelper.DegreesToRadians(-yaw) );
rotation.Transpose(); // OpenTK matrices are transposed by default for some reason
Vector4 rotatedMovement = rotation * new Vector4(movement, 1);
movement.X = -rotatedMovement.X;
movement.Y = rotatedMovement.Y;
movement.Z = -rotatedMovement.Z;
movement *= 0.1f;
foreach (int pID in particleIDs) {
particleSystem.AddVel(pID, movement.X, movement.Y, movement.Z);
}
}float w_poly6(float dist, float h) {
dist = clamp(dist, (float) 0, (float) h);
float tmp = h * h - dist * dist;
return ( 315.0f / (64.0f * PI * pow(h, 9)) ) * tmp * tmp * tmp;
}
float w_spikygrad(float dist, float h) {
dist = clamp(dist, (float) 0, (float) h);
if (dist < FLT_EPSILON) return 0; // too close = same particle = can't use this scalar kernel
return ( -45 / (PI * pow(h, 6)) ) * (h - dist) * (h - dist);
}
kernel void calculate_lambdas(
global float *eposBuffer,
global float *imasses,
global float *lambdas,
global int *cellIDsOfParticles,
global int *cellStartAndEndIDs,
global int *sortedParticleIDs,
global int *phase
) {
int i = get_global_id(0);
if (phase[i] != PHASE_LIQUID) {
lambdas[i] = 0;
return;
}
float3 epos1 = getVec(eposBuffer, i);
float density = 0;
float gradN = 0; // gradient sum when other particle is neighbour
float3 gradS = 0; // gradient sum when other particle is self
FOREACH_NEIGHBOUR_j
float3 epos2 = getVec(eposBuffer, j);
float3 diff = epos1 - epos2;
float dist = length(diff);
if (dist > KERNEL_SIZE) continue;
// the added bit should be multiplied by an extra scalar if its a solid
if (imasses[j] > 0) density += (1.0 / imasses[j]) * w_poly6(dist, KERNEL_SIZE);
if (i != j) {
float kgrad = w_spikygrad(dist, KERNEL_SIZE);
float tmp = kgrad / REST_DENSITY;
// the added bit should be multiplied by an extra scalar if its a solid
gradN += tmp * tmp;
// the added bit should be multiplied by an extra scalar if its a solid
gradS += normalize(diff) * kgrad;
}
END_FOREACH_NEIGHBOUR_j
gradS /= REST_DENSITY;
float denominator = gradN + dot(gradS, gradS);
lambdas[i] = -(density / REST_DENSITY - 1.0) / (denominator + RELAXATION);
}Keybinds:
- W: move forward
- S: move backward
- A: move left
- D: move right
- Space: move up
- Left Shift: move down
- Escape: close program
References: