From cd53dc3fc652be64782016c61898d0ce41cc679d Mon Sep 17 00:00:00 2001
From: Jose Agustin Arguedas Leiva <agustin.arguedas@ds.mpg.de>
Date: Tue, 2 Mar 2021 10:44:47 +0100
Subject: [PATCH 1/2] moved variables into protected to use in derived classes

---
 cpp/particles/p2p/p2p_ghost_collisions.hpp | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/cpp/particles/p2p/p2p_ghost_collisions.hpp b/cpp/particles/p2p/p2p_ghost_collisions.hpp
index 166ac5fc..4671b91a 100644
--- a/cpp/particles/p2p/p2p_ghost_collisions.hpp
+++ b/cpp/particles/p2p/p2p_ghost_collisions.hpp
@@ -45,17 +45,12 @@ class p2p_ghost_collisions
         const double pi  = atan(1.0)*4;
         const double pi2 = atan(1.0)*8;
         // depending on the particle shape, we will be deciding whether particles intersect in different ways
-        enum particle_shape {CYLINDER, SPHERE, DISK};
-        particle_shape current_particle_shape;
-
-        // description for cylinders:
-        double cylinder_width;
-        double cylinder_length;
 
         // description for disks:
         double disk_width;
-        bool isActive;
    protected:        
+       enum particle_shape {CYLINDER, SPHERE, DISK};
+       particle_shape current_particle_shape;
         void add_colliding_pair(partsize_t idx_part1, partsize_t idx_part2)
         {
             // store colliding particle ids in order, to be able to identify pairs more easily
@@ -70,6 +65,12 @@ class p2p_ghost_collisions
             this->collision_pairs_local.push_back(idx_part_small);
             this->collision_pairs_local.push_back(idx_part_big);
         }
+        // description for cylinders:
+        double cylinder_width;
+        double cylinder_length;
+
+        bool isActive;
+        
         bool synchronisation;
         std::vector <partsize_t> collision_pairs_local;
         std::vector <partsize_t> collision_pairs_global;
-- 
GitLab


From 99a1922ed7b765a4c5e8b9c83411b58401ecf00b Mon Sep 17 00:00:00 2001
From: Jose Agustin Arguedas Leiva <agustin.arguedas@ds.mpg.de>
Date: Wed, 3 Mar 2021 13:34:07 +0100
Subject: [PATCH 2/2] passed functionality in p2p cylinders to a function

---
 cpp/particles/p2p/p2p_ghost_collisions.hpp | 55 ++++++++++++----------
 1 file changed, 30 insertions(+), 25 deletions(-)

diff --git a/cpp/particles/p2p/p2p_ghost_collisions.hpp b/cpp/particles/p2p/p2p_ghost_collisions.hpp
index 4671b91a..84bb63b0 100644
--- a/cpp/particles/p2p/p2p_ghost_collisions.hpp
+++ b/cpp/particles/p2p/p2p_ghost_collisions.hpp
@@ -90,6 +90,16 @@ public:
         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>
     void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
     }
@@ -209,7 +219,8 @@ public:
                 break;
             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 y = yseparation;
@@ -218,12 +229,17 @@ public:
 
                     if( r <= this->cylinder_length )
                     {
-
                         /* 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 = 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];
-                        xp = x * pos_part1[IDXC_X+3] + y * pos_part1[IDXC_Y+3] + z * pos_part1[IDXC_Z+3];
-                        xq = x * pos_part2[IDXC_X+3] + y * pos_part2[IDXC_Y+3] + z * pos_part2[IDXC_Z+3];
+                        pq = px * qx + py * qy + pz * qz;
+                        xp = x * px + y * py + z * pz;
+                        xq = x * qx + y * qy + z * qz;
                         /* t and s parametrize the two rods. Find min distance: */
                         assert(this->cylinder_length > 0);
                         t = 2.0/(this->cylinder_length*(pq*pq-1.0))*(xp-pq*xq);
@@ -232,10 +248,7 @@ public:
                         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 = this->cylinder_length*0.5*t*pos_part1[IDXC_X+3]-this->cylinder_length*0.5*s*pos_part2[IDXC_X+3]+x;
-                            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);
+                            min_distance = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
                         }
                         else
                         {
@@ -244,34 +257,26 @@ public:
                             t = 1.0;
                             s = t*pq+2.0/this->cylinder_length*xq;
                             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;
-                            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 = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+                                min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
+                                min_distance = fmin( min_distance_current, min_distance );
                             /* t fixed at -1, find min along s */
                             t = -1.0;
                             s = t*pq+2.0/this->cylinder_length*xq;
                             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;
-                            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 = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+                                min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
+                                min_distance = fmin( min_distance_current, min_distance );
                             /* s fixed at 1, find min along t */
                             s = 1.0;
                             t = s*pq-2.0/this->cylinder_length*xp;
                             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;
-                            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 = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+                                min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
+                                min_distance = fmin( min_distance_current, min_distance );
                             /* s fixed at -1, find min along t */
                             s = -1.0;
                             t = s*pq-2.0/this->cylinder_length*xp;
                             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;
-                            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 = fmin( sqrt(x_dist*x_dist+y_dist*y_dist+z_dist*z_dist), min_distance );
+                                min_distance_current = this->rod_distance(px, py, pz, qx, qy, qz, t, s, x, y, z);
+                                min_distance = fmin( min_distance_current, min_distance );
                         }
                         /* If cylinders overlap count it as a collision */
                         if( min_distance<=this->cylinder_width )
-- 
GitLab