Skip to content
GitLab
Commit 20cc0f06 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

add methods to set p2p_cylinder parameters

parent b1f7b3aa
Hide whitespace changes
Inline Side-by-side
cpp/particles/p2p_cylinder_collisions.hpp
View file @ 20cc0f06
......@@ -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;
}
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment