Skip to content
Snippets Groups Projects
Commit 99a1922e authored by Jose Agustin Arguedas Leiva's avatar Jose Agustin Arguedas Leiva
Browse files

passed functionality in p2p cylinders to a function

parent cd53dc3f
No related branches found
No related tags found
1 merge request!30Feature/write vector
Pipeline #94847 failed
Stage: build
Stage: test
Show whitespace changes
Inline Side-by-side
cpp/particles/p2p/p2p_ghost_collisions.hpp
+ 30
− 25
View file @ 99a1922e
...@@ -90,6 +90,16 @@ public: ...@@ -90,6 +90,16 @@ public:
this->collision_pairs_local.reserve(src.collision_pairs_local.capacity()); this->collision_pairs_local.reserve(src.collision_pairs_local.capacity());
} }
double rod_distance(double px, double py, double pz, double qx, double qy, double qz, double t, double s, double x, double y, double z){
double x_dist, y_dist, z_dist;
double min_distance;
x_dist = this->cylinder_length*0.5*t*px-this->cylinder_length*0.5*s*qx+x;
y_dist = this->cylinder_length*0.5*t*py-this->cylinder_length*0.5*s*qy+y;
z_dist = this->cylinder_length*0.5*t*pz-this->cylinder_length*0.5*s*qz+z;
min_distance = sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist);
return min_distance;
}
template <int size_particle_rhs> template <int size_particle_rhs>
void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{ void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
} }
...@@ -209,7 +219,8 @@ public: ...@@ -209,7 +219,8 @@ public:
break; break;
case CYLINDER: case CYLINDER:
{ {
double pq, xp, xq, t, s, x_dist, y_dist, z_dist ,min_distance; double pq, xp, xq, t, s, min_distance, min_distance_current;
double px,py, pz, qx, qy, qz;
const double x = xseparation; const double x = xseparation;
const double y = yseparation; const double y = yseparation;
...@@ -218,12 +229,17 @@ public: ...@@ -218,12 +229,17 @@ public:
if( r <= this->cylinder_length ) if( r <= this->cylinder_length )
{ {
/* p and q are the orientation vectors of the first and second particles. */ /* p and q are the orientation vectors of the first and second particles. */
px = pos_part1[IDXC_X+3];
py = pos_part1[IDXC_Y+3];
pz = pos_part1[IDXC_Z+3];
qx = pos_part2[IDXC_X+3];
qy = pos_part2[IDXC_Y+3];
qz = pos_part2[IDXC_Z+3];
/* pq, xp, xq are scalar products of these vectors with x, relative position */ /* pq, xp, xq are scalar products of these vectors with x, relative position */
pq = pos_part1[IDXC_X+3] * pos_part2[IDXC_X+3] + pos_part1[IDXC_Y+3] * pos_part2[IDXC_Y+3] + pos_part1[IDXC_Z+3] * pos_part2[IDXC_Z+3]; pq = px * qx + py * qy + pz * qz;
xp = x * pos_part1[IDXC_X+3] + y * pos_part1[IDXC_Y+3] + z * pos_part1[IDXC_Z+3]; xp = x * px + y * py + z * pz;
xq = x * pos_part2[IDXC_X+3] + y * pos_part2[IDXC_Y+3] + z * pos_part2[IDXC_Z+3]; xq = x * qx + y * qy + z * qz;
/* t and s parametrize the two rods. Find min distance: */ /* t and s parametrize the two rods. Find min distance: */
assert(this->cylinder_length > 0); assert(this->cylinder_length > 0);
t = 2.0/(this->cylinder_length*(pq*pq-1.0))*(xp-pq*xq); t = 2.0/(this->cylinder_length*(pq*pq-1.0))*(xp-pq*xq);
...@@ -232,10 +248,7 @@ public: ...@@ -232,10 +248,7 @@ public:
if( abs(t)<=1.0 and abs(s)<=1.0 ) 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 */ /* Get minimal distance in case of both t and s in {-1,1}. Else: check edges */
x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x; min_distance = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y;
z_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Z+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Z+3]+z;
min_distance = sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist);
} }
else else
{ {
...@@ -244,34 +257,26 @@ public: ...@@ -244,34 +257,26 @@ public:
t = 1.0; t = 1.0;
s = t*pq+2.0/this->cylinder_length*xq; s = t*pq+2.0/this->cylinder_length*xq;
if( abs(s)>1.0 ) { s = s / abs(s) ;} if( abs(s)>1.0 ) { s = s / abs(s) ;}
x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x; min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y; min_distance = fmin( min_distance_current, min_distance );
z_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Z+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Z+3]+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 fixed at -1, find min along s */
t = -1.0; t = -1.0;
s = t*pq+2.0/this->cylinder_length*xq; s = t*pq+2.0/this->cylinder_length*xq;
if( abs(s)>1.0 ) { s = s / abs(s) ;} if( abs(s)>1.0 ) { s = s / abs(s) ;}
x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x; min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y; min_distance = fmin( min_distance_current, min_distance );
z_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Z+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Z+3]+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 fixed at 1, find min along t */
s = 1.0; s = 1.0;
t = s*pq-2.0/this->cylinder_length*xp; t = s*pq-2.0/this->cylinder_length*xp;
if( abs(t)>1.0 ) { t = t / abs(t) ;} if( abs(t)>1.0 ) { t = t / abs(t) ;}
x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x; min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y; min_distance = fmin( min_distance_current, min_distance );
z_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Z+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Z+3]+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 fixed at -1, find min along t */
s = -1.0; s = -1.0;
t = s*pq-2.0/this->cylinder_length*xp; t = s*pq-2.0/this->cylinder_length*xp;
if( abs(t)>1.0 ) { t = t / abs(t) ;} if( abs(t)>1.0 ) { t = t / abs(t) ;}
x_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x; min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
y_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Y+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Y+3]+y; min_distance = fmin( min_distance_current, min_distance );
z_dist = this->cylinder_length*0.5*t*pos_part1[IDXC_Z+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_Z+3]+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 cylinders overlap count it as a collision */
if( min_distance<=this->cylinder_width ) if( min_distance<=this->cylinder_width )
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment