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.

Sign up for GitHub

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

Fix bond drawing bugs in the OpenGL visualizer #3511

Merged
merged 4 commits into from
Feb 19, 2020
Merged

Fix bond drawing bugs in the OpenGL visualizer #3511

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
3d6ccab
Fix periodic image bug
jngrad Feb 17, 2020
4d555aa
Remove wildcard import
jngrad Feb 17, 2020
fa36876
Simplify code
jngrad Feb 18, 2020
ec59e22
Fix bond connectivity issues at the box boundaries
jngrad Feb 18, 2020
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
166 changes: 89 additions & 77 deletions src/python/espressomd/visualization_opengl.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
import OpenGL.GLUT
import OpenGL.GLU
import OpenGL.GL
from math import *
import math
import numpy as np
import ctypes
import os
Expand Down Expand Up @@ -413,6 +413,8 @@ class openGLLive:

# HAS PERIODIC IMAGES
self.has_images = any(i != 0 for i in self.specs['periodic_images'])
self.image_vectors = [
list(range(-i, i + 1)) for i in self.specs['periodic_images']]

# INITS
self.system = system
Expand Down Expand Up @@ -1141,15 +1143,13 @@ class openGLLive:
self.particles['pos'][pid], radius, self.specs['quality_particles'])

if self.has_images:
for imx in range(-self.specs['periodic_images'][0],
self.specs['periodic_images'][0] + 1):
for imy in range(-self.specs['periodic_images'][1],
self.specs['periodic_images'][1] + 1):
for imz in range(-self.specs['periodic_images'][2],
self.specs['periodic_images'][2] + 1):
for imx in self.image_vectors[0]:
for imy in self.image_vectors[1]:
for imz in self.image_vectors[2]:
if imx != 0 or imy != 0 or imz != 0:
self._redraw_sphere(
self.particles['pos'][pid] + (imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2]), radius, self.specs['quality_particles'])
offset = [imx, imy, imz] * self.system.box_l
self._redraw_sphere(self.particles['pos'][pid] + offset,
radius, self.specs['quality_particles'])

IF EXTERNAL_FORCES:
if self.specs['ext_force_arrows'] or pid == self.dragId:
Expand Down Expand Up @@ -1195,72 +1195,55 @@ class openGLLive:
v, dtype=float) * sc, radius, col, self.materials['chrome'], self.specs['quality_arrows'])

def _draw_bonds(self):
pIds = range(len(self.particles['pos']))
b2 = self.system.box_l[0] / 2.0
box_l2_sqr = pow(b2, 2.0)
box_l_2 = self.system.box_l / 2.0
for b in self.bonds:
col = self._modulo_indexing(self.specs['bond_type_colors'], b[2])
mat = self.materials[self._modulo_indexing(
self.specs['bond_type_materials'], b[2])]
radius = self._modulo_indexing(
self.specs['bond_type_radius'], b[2])
d = self.particles['pos'][b[0]] - self.particles['pos'][b[1]]
bondLen_sqr = d[0] * d[0] + d[1] * d[1] + d[2] * d[2]
xA = self.particles['pos'][b[0]]
xB = self.particles['pos'][b[1]]
dx = xB - xA

if bondLen_sqr < box_l2_sqr:
if abs(dx[0]) < box_l_2[0] and abs(
dx[1]) < box_l_2[1] and abs(dx[2]) < box_l_2[2]:
# BOND COMPLETELY INSIDE BOX
draw_cylinder(
self.particles['pos'][
b[0]], self.particles['pos'][b[1]], radius,
col, mat, self.specs['quality_bonds'])
for imx in range(-self.specs['periodic_images'][0],
self.specs['periodic_images'][0] + 1):
for imy in range(-self.specs['periodic_images'][1],
self.specs['periodic_images'][1] + 1):
for imz in range(-self.specs['periodic_images'][2],
self.specs['periodic_images'][2] + 1):
if imx != 0 or imy != 0 or imz != 0:
im = np.array([imx, imy, imz])
draw_cylinder(self.particles['pos'][b[0]] + im * self.imPos[dim], self.particles['pos'][b[1]] +
im * self.imPos[dim], radius, col, mat, self.specs['quality_bonds'])
draw_cylinder(xA, xB, radius, col, mat,
self.specs['quality_bonds'])
if self.has_images:
for imx in self.image_vectors[0]:
for imy in self.image_vectors[1]:
for imz in self.image_vectors[2]:
if imx != 0 or imy != 0 or imz != 0:
offset = [imx, imy, imz] * \
self.system.box_l
draw_cylinder(xA + offset, xB + offset,
radius, col, mat, self.specs['quality_bonds'])
else:
# SPLIT BOND
l = self.particles['pos'][b[0]] - self.particles['pos'][b[1]]
l0 = self.particles['pos'][b[0]]
hits = 0
for i in range(6):
lineBoxNDot = float(np.dot(l, self.box_n[i]))
if lineBoxNDot == 0:
continue
s = l0 - \
np.dot(l0 - self.box_p[i],
self.box_n[i]) / lineBoxNDot * l
if self._is_inside_box(s):
if lineBoxNDot < 0:
s0 = s
else:
s1 = s
hits += 1
if hits >= 2:
break
if hits >= 2:
draw_cylinder(self.particles['pos'][b[0]], s0, radius, col,
mat, self.specs['quality_bonds'])
draw_cylinder(self.particles['pos'][b[1]], s1, radius, col,
mat, self.specs['quality_bonds'])

for imx in range(-self.specs['periodic_images'][0],
self.specs['periodic_images'][0] + 1):
for imy in range(-self.specs['periodic_images'][1],
self.specs['periodic_images'][1] + 1):
for imz in range(-self.specs['periodic_images'][2],
self.specs['periodic_images'][2] + 1):
# BOND CROSSES THE BOX BOUNDARIES
d = self._cut_bond(xA, dx)
if d is np.inf:
continue

sA = xA + d * dx
sB = xB - (1 - d) * dx
draw_cylinder(xA, sA, radius, col, mat,
self.specs['quality_bonds'])
draw_cylinder(xB, sB, radius, col, mat,
self.specs['quality_bonds'])

if self.has_images:
for imx in self.image_vectors[0]:
for imy in self.image_vectors[1]:
for imz in self.image_vectors[2]:
if imx != 0 or imy != 0 or imz != 0:
im = np.array([imx, imy, imz])
draw_cylinder(self.particles['pos'][b[0]] + im * self.imPos[dim], s0 + im *
self.imPos[dim], radius, col, mat, self.specs['quality_bonds'])
draw_cylinder(self.particles['pos'][b[1]] + im * self.imPos[dim], s1 + im *
self.imPos[dim], radius, col, mat, self.specs['quality_bonds'])
offset = [imx, imy, imz] * \
self.system.box_l
draw_cylinder(xA + offset, sA + offset, radius, col, mat,
self.specs['quality_bonds'])
draw_cylinder(xB + offset, sB + offset, radius, col, mat,
self.specs['quality_bonds'])

def _redraw_sphere(self, pos, radius, quality):
OpenGL.GL.glPushMatrix()
Expand All @@ -1269,12 +1252,44 @@ class openGLLive:
OpenGL.GL.OpenGL.GL.glPopMatrix()

# HELPER TO DRAW PERIODIC BONDS
def _is_inside_box(self, p):
eps = 1e-5
def _cut_bond(self, xA, dx):
"""
Algorithm for the line-plane intersection problem. Given the unfolded
positions of two particles, determine: 1) the box image that minimizes
the folded bond length, 2) which side of the simulation box is crossed
by the bond, 3) how much of the bond is inside the unit cell starting
from the first particle. That scalar is returned by the function, and
can be used to calculate the coordinates of the line-plane
intersection point. The algorithm can be found at
https://en.wikipedia.org/w/index.php?title=Line%E2%80%93plane_intersection&oldid=940427752#Algebraic_form
"""
# corner case: the two particles occupy the same position
if np.dot(dx, dx) < 1e-9:
return np.inf
# find the box image that minimizes the unfolded bond length
shift = np.rint(dx / self.system.box_l)
# get the unfolded bond vector
dx -= shift * self.system.box_l
# find which side of the simulation box is crossed by the bond
best_d = np.inf
best_i = 0
for i in range(3):
if p[i] < -eps or p[i] > eps + self.last_box_l[i]:
return False
return True
if dx[i] == 0:
continue # corner case: bond is parallel to a face
elif dx[i] > 0:
p0_i = self.system.box_l[i]
else:
p0_i = 0
# calculate the length of the bond that is inside the box using
# an optimized version of `np.dot(p0 - x0, n) / np.dot(dx, n)`
# where the dot products decay to an array access, because `n`
# is always a unit vector in rectangular boxes and its sign
# is canceled out by the division
d = (p0_i - xA[i]) / dx[i]
if d < best_d:
best_d = d
best_i = i
return best_d

# ARROWS IN A PLANE FOR LB VELOCITIES
def _draw_lb_vel(self):
Expand Down Expand Up @@ -1619,9 +1634,6 @@ class openGLLive:

self.depth = 0

self.imPos = [np.array([self.system.box_l[0], 0, 0]), np.array(
[0, self.system.box_l[1], 0]), np.array([0, 0, self.system.box_l[2]])]

# LOOK FOR LB ACTOR
if self.specs['LB_draw_velocity_plane'] or self.specs['LB_draw_nodes'] or self.specs['LB_draw_node_boundaries']:
lb_types = [espressomd.lb.LBFluid]
Expand Down Expand Up @@ -2021,7 +2033,7 @@ def rotation_helper(d):
# the rotation axis is the cross product between z and d
vz = np.cross([0.0, 0.0, 1.0], d)
# get the angle using a dot product
angle = 180.0 / np.pi * acos(d[2] / np.linalg.norm(d))
angle = 180.0 / np.pi * math.acos(d[2] / np.linalg.norm(d))

return angle, vz[0], vz[1]

Expand Down Expand Up @@ -2218,7 +2230,7 @@ def draw_sphero_cylinder(posA, posB, radius, color, material, quality):
if v == 0:
ax = 57.2957795
else:
ax = 57.2957795 * acos(d[2] / v)
ax = 57.2957795 * math.acos(d[2] / v)

if d[2] < 0.0:
ax = -ax
Expand Down Expand Up @@ -2575,8 +2587,8 @@ class Camera:
return m

def rotate_vector(self, vec, ang, axe):
sinhalf = sin(ang / 2)
coshalf = cos(ang / 2)
sinhalf = math.sin(ang / 2)
coshalf = math.cos(ang / 2)

rx = axe[0] * sinhalf
ry = axe[1] * sinhalf
Expand Down