Update ExternalFieldWall.h

Added cylinder-spheropolyhedra overlap method

Author: seoulfood

hoomd/hpmc/ExternalFieldWall.h

@@ -481,6 +481,86 @@ inline bool test_confined(const CylinderWall& wall,
     return accept;
     }
 
+// Cylindrical Walls and Convex Spheropolyhedra
+DEVICE inline bool test_confined(const CylinderWall& wall,
+                                 const ShapeSpheropolyhedron& shape,
+                                 const vec3<Scalar>& position,
+                                 const vec3<Scalar>& box_origin,
+                                 const BoxDim& box)
+    {
+    bool accept = true;
+    Scalar3 t = vec_to_scalar3(position - box_origin);
+    t.x = t.x - wall.origin.x;
+    t.y = t.y - wall.origin.y;
+    t.z = t.z - wall.origin.z;
+    vec3<Scalar> shifted_pos(box.minImage(t));
+
+    vec3<Scalar> dist_vec
+        = cross(shifted_pos,
+                wall.orientation); // find the component of the shifted position that is
+                                   // perpendicular to the normalized orientation vector
+    Scalar max_dist = sqrt(dot(dist_vec, dist_vec));
+    if (wall.inside)
+        {
+        max_dist += (shape.getCircumsphereDiameter() / Scalar(2.0));        
+        }
+    else
+        {
+        max_dist -= (shape.getCircumsphereDiameter() / Scalar(2.0));
+        if (max_dist < 0)
+            {
+            max_dist = 0;
+            }
+
+        }
+
+    bool check_verts
+        = wall.inside ? (wall.rsq <= max_dist * max_dist)
+                      : (wall.rsq >= max_dist * max_dist); // condition to check vertices, dependent
+                                                           // on inside or outside container
+
+    if (check_verts)
+        {
+        if (shape.verts.N)
+            {
+            if (wall.inside)
+                {
+                for (unsigned int v = 0; v < shape.verts.N && accept; v++)
+                    {
+                    vec3<Scalar> pos(shape.verts.x[v], shape.verts.y[v], shape.verts.z[v]);
+                    vec3<Scalar> rotated_pos = rotate(shape.orientation, pos);
+                    rotated_pos += shifted_pos;
+                    
+                    dist_vec = cross(rotated_pos, wall.orientation);
+                    max_dist = sqrt(dot(dist_vec, dist_vec));
+                    Scalar max_tot_dist = max_dist + shape.verts.sweep_radius;
+
+                    accept = wall.inside ? (wall.rsq > max_tot_dist * max_tot_dist)
+                                        : (wall.rsq < max_tot_dist * max_tot_dist);
+                    }
+                }
+            else
+                {
+                // build a sphero-polyhedron and for the wall and the convex polyhedron
+                quat<Scalar> q; // default is (1, 0, 0, 0)
+                unsigned int err = 0;
+                ShapeSpheropolyhedron wall_shape(q, *wall.verts);
+                ShapeSpheropolyhedron part_shape(shape.orientation, shape.verts);
+
+                accept = !test_overlap(shifted_pos, wall_shape, part_shape, err);
+                }
+            }
+        // Edge case; pure sphere. In this case, check_verts implies that the
+        // sphere will be outside. 
+        // Note from gabs: I don't fully understand this edge case
+        else
+            {
+            accept = false;
+            }
+        }
+    return accept;
+    }
+
 // Plane Walls and Spheres
 template<>
 inline bool test_confined<PlaneWall, ShapeSphere>(const PlaneWall& wall,