Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Cylinder-Wall Spheropolyhedron Support #1984

Open
wants to merge 3 commits into
base: trunk-patch
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
80 changes: 80 additions & 0 deletions hoomd/hpmc/ExternalFieldWall.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Comment on lines +498 to +501
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Comments are easier to read at the same level as the code. You can also use English language rules when writing comments.

Suggested change
vec3<Scalar> dist_vec
= cross(shifted_pos,
wall.orientation); // find the component of the shifted position that is
// perpendicular to the normalized orientation vector
// Find the component of the shifted position that is
// perpendicular to the normalized orientation vector.
vec3<Scalar> dist_vec
= cross(shifted_pos,
wall.orientation);

Scalar max_dist = sqrt(dot(dist_vec, dist_vec));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This value is only ever used squared.

Remove the sqrt and rename to max_dist_squared. Remove the multiplications that square it.

if (wall.inside)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A comment would be helpful to understand what the code is doing here.

Suggested change
if (wall.inside)
// Check if the shape's bounding sphere intersects the wall.
if (wall.inside)

Also, consider rewriting this so that it evaluates the <= or > inside the if. The branch is repeated again with a ternary operator on the next line which is very confusing. Better yet, combine the content of the if check_verts -> if shape.verts.N -> branches here as well. Bringing all the "inside" and "outside" code together will make it much easier to understand.

{
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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When shape.verts.N is zero, it is assume to have a single vertex at 0,0,0 - making it a sphere. In that case the bounding sphere check above is the only check that is necessary. I agree that this is really confusing. This case should be handled FIRST, not last. It would be more readable moved up above and handled with its own return false branch.

else
{
accept = false;
}
}
return accept;
}

// Plane Walls and Spheres
template<>
inline bool test_confined<PlaneWall, ShapeSphere>(const PlaneWall& wall,
Expand Down
5 changes: 4 additions & 1 deletion hoomd/hpmc/external/wall.py
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,10 @@ class _HPMCWallsMetaList(_WallsMetaList):
hoomd.wall.Cylinder,
hoomd.wall.Plane,
],
hpmc.integrate.ConvexSpheropolyhedron: [hoomd.wall.Sphere, hoomd.wall.Plane],
hpmc.integrate.ConvexSpheropolyhedron: [
hoomd.wall.Sphere,
hoomd.wall.Cylinder,
hoomd.wall.Plane],
}

def _check_wall_compatibility(self, wall):
Expand Down
3 changes: 1 addition & 2 deletions hoomd/hpmc/integrate.py
Original file line number Diff line number Diff line change
Expand Up @@ -1676,8 +1676,7 @@ class ConvexSpheropolyhedron(HPMCIntegrator):
.. rubric:: Wall support.
`ConvexSpheropolyhedron` supports the `hoomd.wall.Sphere` and
`hoomd.wall.Plane` geometries.
`ConvexSpheropolyhedron` supports all `hoomd.wall` geometries.
Example::
Expand Down
Loading