I am trying to implement a particle-AABOX edge collision, the below images represent, two timesteps(dt) where the spherical particle is accelerated by gravity.
I have the particle center C and its radius r, along with AABOX coordinates as B_min and Bmax.
I have to check if the particle collides with any edges, and then bounce back, if a collision happens.
If I use Arvo's AABB vs Sphere collision algorithm, then I can find if there is an intersection within the radius, but I can't get the hitpoint and hitdistance.
And after getting the hitpoint and hitdistance, how to reflect it back? (for every timestep, I will have the current position C, and also the previous position P).
I tried to calculate it by implementing,
BOX_N[0] = make_vector(-1, 0, 0);// x-min;
BOX_P[0] = make_point(bb_min.x, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_min.z + (bb_max.z - bb_min.z) / 2.0f);
BOX_N[1] = make_vector(1, 0, 0); // x-max
BOX_P[1] = make_point(bb_max.x, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_min.z + (bb_max.z - bb_min.z) / 2.0f);
BOX_N[2] = make_vector(0, -1, 0);// y-min;
BOX_P[2] = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_min.y, bb_min.z + (bb_max.z - bb_min.z) / 2.0f);
BOX_N[3] = make_vector(0, 1, 0); // y-max
BOX_P[3] = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_max.y, bb_min.z + (bb_max.z - bb_min.z) / 2.0f);
BOX_N[4] = make_vector(0, 0, -1);// z-min;
BOX_P[4] = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_min.z);
BOX_N[5] = make_vector(0, 0, 1); // z-max
BOX_P[5] = make_point(bb_min.x + (bb_max.x - bb_min.x) / 2.0f, bb_min.y + (bb_max.y - bb_min.y) / 2.0f , bb_max.z);
__device__ int box_collision(const Point& previous_position, const Point& current_position, const Vector& direction, const float& radius, Point& hitPoint, Vector& normal, float& hit_distance) {
Point boxPoint;
int index;
if (current_position.x < (bb_min.x + radius)) { index = 1; normal = BOX_N[0]; boxPoint = BOX_P[0]; boxPoint.x += radius;}
else if (current_position.x > (bb_max.x - radius)) { index = 1; normal = BOX_N[1]; boxPoint = BOX_P[1]; boxPoint.x -= radius;}
else if (current_position.y < (bb_min.y + radius)) { index = 2; normal = BOX_N[2]; boxPoint = BOX_P[2]; boxPoint.y += radius;}
else if (current_position.y > (bb_max.y - radius)) { index = 2; normal = BOX_N[3]; boxPoint = BOX_P[3]; boxPoint.y -= radius;}
else if (current_position.z < (bb_min.z + radius)) { index = 3; normal = BOX_N[4]; boxPoint = BOX_P[4]; boxPoint.z += radius;}
else if (current_position.z > (bb_max.z - radius)) { index = 3; normal = BOX_N[5]; boxPoint = BOX_P[5]; boxPoint.z -= radius;}
else return 0;
auto denom = vdot(direction, normal);
if (denom < 1e-6) return 0;
hit_distance = vdot(normal, boxPoint - previous_position) / denom;
if (hit_distance < 0) return 0;
hitPoint = previous_position + hit_distance * direction;
return index;
}
inline __device__ bool compute_box_collision(Point& previous_position, Point& current_position, const float& radius) {
Point hitPoint;
Vector normal;
float hit_distance;
Vector direction = vnormalize(current_position - previous_position);
int collision = box_collision(previous_position, current_position, direction, radius, hitPoint, normal, hit_distance);
if (!collision) return true;
Vector damping{1, 1, 1};
if (collision == 1) damping.x = bounce_factor;
else if (collision == 2) damping.y = bounce_factor;
else if (collision == 3) damping.z = bounce_factor;
float relection_length = vlength(hitPoint - current_position);
Vector R = vnormalize(direction - 2.0f * vdot(normal, direction) * normal) * damping;
current_position = hitPoint + R * relection_length;
previous_position = hitPoint - R * hit_distance;
return false;
}
by calling the 2nd function to check for collision,
inline __device__ bool compute_box_collision(Point& previous_position, Point& current_position, const float& radius)
In the above code, I used the ray-plane interesection, and then use reflection of the direction to get the bounce direction, but after 500-550 timesteps, the solution deviates and becomes numerically unstable(possibly they settle at the bottom and collides continously). Is there something, I am missing here? Or am I doing something wrong?
Note: I am using simiulating 10k particles, so the solution deviation compounds by a large amount.