diff --git a/cpp/particles/p2p_cylinder_collisions.hpp b/cpp/particles/p2p_cylinder_collisions.hpp index dc4bc735899e13b66262a5990c45a60f7f4b0ffa..f68b4d3ec8e5075ca2214d5b31e7035b1c47c013 100644 --- a/cpp/particles/p2p_cylinder_collisions.hpp +++ b/cpp/particles/p2p_cylinder_collisions.hpp @@ -32,10 +32,10 @@ class p2p_cylinder_ghost_collisions{ double width; double length; /* Following doubles are needed for the collision computation */ - + public: - p2p_cylinder_ghost_collisions() : collision_counter(0){} + p2p_cylinder_ghost_collisions() : collision_counter(0),width(1.),length(1.){} // Copy constructor use a counter set to zero p2p_cylinder_ghost_collisions(const p2p_cylinder_ghost_collisions&) : collision_counter(0){} @@ -55,60 +55,60 @@ public: const real_number pos_part2[], real_number /*rhs_part2*/[], const real_number /*dist_pow2*/, const real_number /*cutoff*/, const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/){ - double x, y, z, pq, xp, xq, t, s, x_dist, y_dist, z_dist ,min_distance; - - /* Relative position vector of the two particles (x,y,z)^T */ - x = pos_part1[0]-pos_part2[0]; - y = pos_part1[1]-pos_part2[1]; - z = pos_part1[2]-pos_part2[2]; - /* p and q are the orientation vectors of the first and second particles. */ - /* pq, xp, xq are scalar products of these vectors with x, relative position */ - pq = pos_part1[3] * pos_part2[3] + pos_part1[4] * pos_part2[4] + pos_part1[5] * pos_part2[5]; - xp = x * pos_part1[3] + y * pos_part1[4] + z * pos_part1[5]; - xq = x * pos_part2[3] + y * pos_part2[4] + z * pos_part2[5]; - /* t and s parametrize the two rods. Find min distance: */ - t = 2.0/(length*(pq*pq-1.0))*(-xp+pq*xq); - s = 2.0/(length*(pq*pq-1.0))*(-pq*xp+xq); - /* Test if -1<s<1 and -1<t<1 */ - if( abs(t)<=1.0 and abs(s)<=1.0 ) { - /* Get minimal distance in case of both t and s in {-1,1}. Else: check edges */ - x_dist = length*0.5*t*pos_part1[3]-length*0.5*s*pos_part2[3]-x; - y_dist = length*0.5*t*pos_part1[4]-length*0.5*s*pos_part2[4]-y; - z_dist = length*0.5*t*pos_part1[5]-length*0.5*s*pos_part2[5]-z; - min_distance = sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist); - } else { - min_distance = length; - /* t fixed at 1, find min along s */ - t = 1.0; - s = t*pq-2.0/length*xq; - x_dist = length*0.5*t*pos_part1[3]-length*0.5*s*pos_part2[3]-x; - y_dist = length*0.5*t*pos_part1[4]-length*0.5*s*pos_part2[4]-y; - z_dist = length*0.5*t*pos_part1[5]-length*0.5*s*pos_part2[5]-z; - min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance ); - /* t fixed at -1, find min along s */ - t = -1.0; - s = t*pq-2.0/length*xq; - x_dist = length*0.5*t*pos_part1[3]-length*0.5*s*pos_part2[3]-x; - y_dist = length*0.5*t*pos_part1[4]-length*0.5*s*pos_part2[4]-y; - z_dist = length*0.5*t*pos_part1[5]-length*0.5*s*pos_part2[5]-z; - min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance ); - /* s fixed at 1, find min along t */ - s = 1.0; - t = s*pq+2.0/length*xp; - x_dist = length*0.5*t*pos_part1[3]-length*0.5*s*pos_part2[3]-x; - y_dist = length*0.5*t*pos_part1[4]-length*0.5*s*pos_part2[4]-y; - z_dist = length*0.5*t*pos_part1[5]-length*0.5*s*pos_part2[5]-z; - min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance ); - /* s fixed at -1, find min along t */ - s = -1.0; - t = s*pq+2.0/length*xp; - x_dist = length*0.5*t*pos_part1[3]-length*0.5*s*pos_part2[3]-x; - y_dist = length*0.5*t*pos_part1[4]-length*0.5*s*pos_part2[4]-y; - z_dist = length*0.5*t*pos_part1[5]-length*0.5*s*pos_part2[5]-z; - min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance ); - } - /* If cylinders overlap count it as a collision */ - if( min_distance<=width ) + double x, y, z, pq, xp, xq, t, s, x_dist, y_dist, z_dist ,min_distance; + + /* Relative position vector of the two particles (x,y,z)^T */ + x = pos_part1[0]-pos_part2[0]; + y = pos_part1[1]-pos_part2[1]; + z = pos_part1[2]-pos_part2[2]; + /* p and q are the orientation vectors of the first and second particles. */ + /* pq, xp, xq are scalar products of these vectors with x, relative position */ + pq = pos_part1[3] * pos_part2[3] + pos_part1[4] * pos_part2[4] + pos_part1[5] * pos_part2[5]; + xp = x * pos_part1[3] + y * pos_part1[4] + z * pos_part1[5]; + xq = x * pos_part2[3] + y * pos_part2[4] + z * pos_part2[5]; + /* t and s parametrize the two rods. Find min distance: */ + t = 2.0/(length*(pq*pq-1.0))*(-xp+pq*xq); + s = 2.0/(length*(pq*pq-1.0))*(-pq*xp+xq); + /* Test if -1<s<1 and -1<t<1 */ + if( abs(t)<=1.0 and abs(s)<=1.0 ) { + /* Get minimal distance in case of both t and s in {-1,1}. Else: check edges */ + x_dist = length*0.5*t*pos_part1[3]-length*0.5*s*pos_part2[3]-x; + y_dist = length*0.5*t*pos_part1[4]-length*0.5*s*pos_part2[4]-y; + z_dist = length*0.5*t*pos_part1[5]-length*0.5*s*pos_part2[5]-z; + min_distance = sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist); + } else { + min_distance = length; + /* t fixed at 1, find min along s */ + t = 1.0; + s = t*pq-2.0/length*xq; + x_dist = length*0.5*t*pos_part1[3]-length*0.5*s*pos_part2[3]-x; + y_dist = length*0.5*t*pos_part1[4]-length*0.5*s*pos_part2[4]-y; + z_dist = length*0.5*t*pos_part1[5]-length*0.5*s*pos_part2[5]-z; + min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance ); + /* t fixed at -1, find min along s */ + t = -1.0; + s = t*pq-2.0/length*xq; + x_dist = length*0.5*t*pos_part1[3]-length*0.5*s*pos_part2[3]-x; + y_dist = length*0.5*t*pos_part1[4]-length*0.5*s*pos_part2[4]-y; + z_dist = length*0.5*t*pos_part1[5]-length*0.5*s*pos_part2[5]-z; + min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance ); + /* s fixed at 1, find min along t */ + s = 1.0; + t = s*pq+2.0/length*xp; + x_dist = length*0.5*t*pos_part1[3]-length*0.5*s*pos_part2[3]-x; + y_dist = length*0.5*t*pos_part1[4]-length*0.5*s*pos_part2[4]-y; + z_dist = length*0.5*t*pos_part1[5]-length*0.5*s*pos_part2[5]-z; + min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance ); + /* s fixed at -1, find min along t */ + s = -1.0; + t = s*pq+2.0/length*xp; + x_dist = length*0.5*t*pos_part1[3]-length*0.5*s*pos_part2[3]-x; + y_dist = length*0.5*t*pos_part1[4]-length*0.5*s*pos_part2[4]-y; + z_dist = length*0.5*t*pos_part1[5]-length*0.5*s*pos_part2[5]-z; + min_distance = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance ); + } + /* If cylinders overlap count it as a collision */ + if( min_distance<=width ) collision_counter += 1; } @@ -127,6 +127,26 @@ public: void reset_collision_counter(){ collision_counter = 0; } + + void set_width(const double WIDTH) + { + this->width = WIDTH; + } + + double get_width() + { + return this->width; + } + + void set_length(const double LENGTH) + { + this->length = LENGTH; + } + + double get_length() + { + return this->length; + } };