diff --git a/doc/sphinx/constraints.rst b/doc/sphinx/constraints.rst index 4736fd6da82..115c4826f1a 100644 --- a/doc/sphinx/constraints.rst +++ b/doc/sphinx/constraints.rst @@ -50,6 +50,7 @@ Available shapes are listed below. - :class:`espressomd.shapes.Sphere` - :class:`espressomd.shapes.SpheroCylinder` - :class:`espressomd.shapes.Stomatocyte` + - :class:`espressomd.shapes.Torus` - :class:`espressomd.shapes.HollowConicalFrustum` - :class:`espressomd.shapes.Union` @@ -400,13 +401,12 @@ The region is described as a pore (lower vertical part of the "T"-shape) and a c .. figure:: figures/slitpore.png :alt: Schematic for the Slitpore shape with labeled geometrical parameters. :align: center - :height: 6.00000cm + :height: 10.00000cm The parameter ``channel_width`` specifies the distance between the top and the plateau edge. The parameter ``pore_length`` specifies the distance between the bottom and the plateau edge. The parameter ``pore_width`` specifies the distance between the two plateau edges, it is the space between the left and right walls of the pore region. -The parameter ``pore_mouth`` specifies the location (z-coordinate) of the pore opening (center). It is always centered in the x-direction. -The parameter ``dividing_plane`` specifies the location (z-coordinate) of the middle between the two walls. +The parameters ``pore_mouth`` and ``dividing_plane`` specify the location in the z-coordinate resp. x-coordinate of the pore opening. All the edges are smoothed via the parameters ``upper_smoothing_radius`` (for the concave corner at the edge of the plateau region) and ``lower_smoothing_radius`` (for the convex corner at the bottom of the pore region). The meaning of the geometrical parameters can be inferred from the schematic in Fig. :ref:`slitpore `. @@ -427,7 +427,7 @@ Pictured is an example constraint with a ``Slitpore`` shape created with :: pore_length=20, pore_mouth=30, pore_width=5, - dividing_plane=40) + dividing_plane=25) system.constraints.add(shape=slitpore, particle_type=0, penetrable=True) @@ -459,6 +459,16 @@ Pictured is an example constraint with a ``SpheroCylinder`` shape created with : system.constraints.add(shape=spherocylinder, particle_type=0) +Torus +""""" + +:class:`espressomd.shapes.Torus` + +It is positioned at ``center`` and has a radius ``radius`` with tube radius ``tube_radius``. +The ``normal`` parameter is the torus rotation axis, which is normalized in the program. +The direction ``direction`` determines the force direction, ``-1`` for inward and ``+1`` for outward. + + HollowConicalFrustum """""""""""""""""""" diff --git a/doc/sphinx/figures/slitpore.pdf b/doc/sphinx/figures/slitpore.pdf deleted file mode 100644 index 874d751fa71..00000000000 Binary files a/doc/sphinx/figures/slitpore.pdf and /dev/null differ diff --git a/doc/sphinx/figures/slitpore.png b/doc/sphinx/figures/slitpore.png index 409d2093912..2b156ac84f5 100644 Binary files a/doc/sphinx/figures/slitpore.png and b/doc/sphinx/figures/slitpore.png differ diff --git a/doc/sphinx/figures/slitpore.svg b/doc/sphinx/figures/slitpore.svg new file mode 100644 index 00000000000..d41f3877abd --- /dev/null +++ b/doc/sphinx/figures/slitpore.svg @@ -0,0 +1,521 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + Copyright (C) 2020 The ESPResSo project + + + + + + + channel_width + + z-axis + + + + + + + + pore_width + lower_smoothing_radius + upper_smoothing_radius + x-axis + + (0, 0) + + + dividing_plane + pore_mouth + + pore_length + + + + diff --git a/samples/visualization_constraints.py b/samples/visualization_constraints.py index 82683480a63..28b092d50c5 100644 --- a/samples/visualization_constraints.py +++ b/samples/visualization_constraints.py @@ -29,7 +29,7 @@ group = parser.add_mutually_exclusive_group() group.add_argument("--wall", action="store_const", dest="shape", const="Wall", default="Wall") -for shape in ("Sphere", "Ellipsoid", "Cylinder", "SpheroCylinder", +for shape in ("Sphere", "Ellipsoid", "Cylinder", "SpheroCylinder", "Torus", "Stomatocyte", "SimplePore", "Slitpore", "HollowConicalFrustum"): group.add_argument("--" + shape.lower(), action="store_const", dest="shape", const=shape) @@ -60,7 +60,7 @@ if args.shape == "Wall": system.constraints.add(shape=espressomd.shapes.Wall( - dist=20, normal=[0.1, 0.0, 1]), + dist=20, normal=[0, 0, 1]), particle_type=0, penetrable=True) elif args.shape == "Sphere": @@ -98,9 +98,9 @@ elif args.shape == "Slitpore": system.constraints.add(shape=espressomd.shapes.Slitpore( - channel_width=15, lower_smoothing_radius=3, upper_smoothing_radius=3, - pore_length=20, pore_mouth=30, pore_width=5), particle_type=0, - penetrable=True) + channel_width=15, lower_smoothing_radius=2, upper_smoothing_radius=3, + pore_length=20, pore_mouth=30, pore_width=10, dividing_plane=25), + particle_type=0, penetrable=True) elif args.shape == "HollowConicalFrustum": system.constraints.add(shape=espressomd.shapes.HollowConicalFrustum( @@ -108,11 +108,22 @@ axis=[0.0, 1.0, 1.0], center=[25, 25, 25], direction=1), particle_type=0, penetrable=True) +elif args.shape == "Torus": + system.constraints.add( + shape=espressomd.shapes.Torus(center=[25] * 3, normal=[1, 1, 1], + direction=1, radius=15, tube_radius=6), + particle_type=0, penetrable=True) + else: raise ValueError("Unknown shape '{}'".format(args.shape)) for i in range(100): - rpos = np.random.random(3) * box_l + # place particles outside the shape + while True: + rpos = np.random.random(3) * box_l + dist = system.constraints[0].shape.calc_distance(position=rpos) + if dist[0] > 2.5: + break system.part.add(pos=rpos, type=1) system.non_bonded_inter[1, 1].lennard_jones.set_params( @@ -126,7 +137,7 @@ max_displacement=0.05) system.integrator.run(500) system.integrator.set_vv() -system.thermostat.set_langevin(kT=10.0, gamma=10, seed=42) +system.thermostat.set_langevin(kT=10.0, gamma=0.1, seed=42) system.force_cap = 1000.0 diff --git a/src/python/espressomd/shapes.py b/src/python/espressomd/shapes.py index a6baceb033b..2d2c5fa2394 100644 --- a/src/python/espressomd/shapes.py +++ b/src/python/espressomd/shapes.py @@ -100,7 +100,10 @@ class Rhomboid(Shape, ScriptInterfaceHelper): class Slitpore(Shape, ScriptInterfaceHelper): """ - .. image:: figures/slitpore.png + .. figure:: figures/slitpore.png + :alt: Schematic for the Slitpore shape with labeled geometrical parameters. + :align: center + :height: 8.00000cm Attributes ---------- diff --git a/src/python/espressomd/visualization_opengl.py b/src/python/espressomd/visualization_opengl.py index d454c14c1b3..41ae100be8b 100644 --- a/src/python/espressomd/visualization_opengl.py +++ b/src/python/espressomd/visualization_opengl.py @@ -2139,6 +2139,8 @@ def __init__(self, shape, particle_type, color, material, self.pore_length = np.array(self.shape.get_parameter('pore_length')) self.pore_mouth = np.array(self.shape.get_parameter('pore_mouth')) self.pore_width = np.array(self.shape.get_parameter('pore_width')) + self.dividing_plane = np.array( + self.shape.get_parameter('dividing_plane')) self.max_box_l = max(self.box_l) def draw(self): @@ -2146,13 +2148,14 @@ def draw(self): # If pore is large, an additional wall is necessary if self.pore_width > 2. * self.lower_smoothing_radius: wall_0 = [ - [0.5 * (self.max_box_l - self.pore_width) + self.lower_smoothing_radius, + [self.dividing_plane - 0.5 * self.pore_width + self.lower_smoothing_radius, 0., self.pore_mouth - self.pore_length], - [0.5 * (self.max_box_l + self.pore_width) - self.lower_smoothing_radius, + [self.dividing_plane + 0.5 * self.pore_width - self.lower_smoothing_radius, 0., self.pore_mouth - self.pore_length], - [0.5 * (self.max_box_l + self.pore_width) - self.lower_smoothing_radius, + [self.dividing_plane + 0.5 * self.pore_width - self.lower_smoothing_radius, self.max_box_l, self.pore_mouth - self.pore_length], - [0.5 * (self.max_box_l - self.pore_width) + self.lower_smoothing_radius, self.max_box_l, self.pore_mouth - self.pore_length]] + [self.dividing_plane - 0.5 * self.pore_width + self.lower_smoothing_radius, + self.max_box_l, self.pore_mouth - self.pore_length]] draw_plane(wall_0, self.color, self.material) # Add the remaining walls @@ -2164,36 +2167,39 @@ def draw(self): wall_2 = [ [0., 0., self.pore_mouth], - [0.5 * (self.max_box_l - self.pore_width) - + [self.dividing_plane - 0.5 * self.pore_width - self.upper_smoothing_radius, 0., self.pore_mouth], - [0.5 * (self.max_box_l - self.pore_width) - + [self.dividing_plane - 0.5 * self.pore_width - self.upper_smoothing_radius, self.max_box_l, self.pore_mouth], [0., self.max_box_l, self.pore_mouth]] wall_3 = [ - [0.5 * (self.max_box_l + self.pore_width) + + [self.dividing_plane + 0.5 * self.pore_width + self.upper_smoothing_radius, 0., self.pore_mouth], [self.max_box_l, 0., self.pore_mouth], [self.max_box_l, self.max_box_l, self.pore_mouth], - [0.5 * (self.max_box_l + self.pore_width) + self.upper_smoothing_radius, self.max_box_l, self.pore_mouth]] + [self.dividing_plane + 0.5 * self.pore_width + + self.upper_smoothing_radius, self.max_box_l, self.pore_mouth]] wall_4 = [ - [0.5 * (self.max_box_l - self.pore_width), 0., + [self.dividing_plane - 0.5 * self.pore_width, 0., self.pore_mouth - self.upper_smoothing_radius], - [0.5 * (self.max_box_l - self.pore_width), self.max_box_l, + [self.dividing_plane - 0.5 * self.pore_width, self.max_box_l, self.pore_mouth - self.upper_smoothing_radius], - [0.5 * (self.max_box_l - self.pore_width), self.max_box_l, self.pore_mouth - - self.pore_length + self.lower_smoothing_radius], - [0.5 * (self.max_box_l - self.pore_width), 0., self.pore_mouth - self.pore_length + self.lower_smoothing_radius]] + [self.dividing_plane - 0.5 * self.pore_width, self.max_box_l, + self.pore_mouth - self.pore_length + self.lower_smoothing_radius], + [self.dividing_plane - 0.5 * self.pore_width, 0., + self.pore_mouth - self.pore_length + self.lower_smoothing_radius]] wall_5 = [ - [0.5 * (self.max_box_l + self.pore_width), 0., + [self.dividing_plane + 0.5 * self.pore_width, 0., self.pore_mouth - self.upper_smoothing_radius], - [0.5 * (self.max_box_l + self.pore_width), self.max_box_l, + [self.dividing_plane + 0.5 * self.pore_width, self.max_box_l, self.pore_mouth - self.upper_smoothing_radius], - [0.5 * (self.max_box_l + self.pore_width), self.max_box_l, self.pore_mouth - - self.pore_length + self.lower_smoothing_radius], - [0.5 * (self.max_box_l + self.pore_width), 0., self.pore_mouth - self.pore_length + self.lower_smoothing_radius]] + [self.dividing_plane + 0.5 * self.pore_width, self.max_box_l, + self.pore_mouth - self.pore_length + self.lower_smoothing_radius], + [self.dividing_plane + 0.5 * self.pore_width, 0., + self.pore_mouth - self.pore_length + self.lower_smoothing_radius]] draw_plane(wall_1, self.color, self.material) draw_plane(wall_2, self.color, self.material) @@ -2206,7 +2212,7 @@ def draw(self): OpenGL.GL.glPushMatrix() quadric = OpenGL.GLU.gluNewQuadric() - OpenGL.GL.glTranslate(0.5 * self.max_box_l - self.upper_smoothing_radius - + OpenGL.GL.glTranslate(self.dividing_plane - self.upper_smoothing_radius - 0.5 * self.pore_width, 0, self.pore_mouth - self.upper_smoothing_radius) OpenGL.GL.glRotatef(ax, rx, ry, 0.) diff --git a/src/shapes/src/Torus.cpp b/src/shapes/src/Torus.cpp index b0416cd7d12..634ad505865 100644 --- a/src/shapes/src/Torus.cpp +++ b/src/shapes/src/Torus.cpp @@ -31,9 +31,9 @@ void Torus::calculate_dist(const Utils::Vector3d &pos, double &dist, auto const r_vec = c_dist - z * e_z; auto const r = r_vec.norm(); - dist = (sqrt(Utils::sqr(r - m_rad) + z * z) - m_tube_rad) * m_direction; + dist = (std::sqrt(Utils::sqr(r - m_rad) + z * z) - m_tube_rad) * m_direction; Utils::Vector3d const dir_vec = c_dist - r_vec * m_rad / r; auto const dir_vec_norm = dir_vec / dir_vec.norm(); - vec = -dir_vec_norm * dist; + vec = dir_vec_norm * std::abs(dist); } } // namespace Shapes diff --git a/testsuite/python/constraint_shape_based.py b/testsuite/python/constraint_shape_based.py index 242f96cf680..a58b9fbf610 100644 --- a/testsuite/python/constraint_shape_based.py +++ b/testsuite/python/constraint_shape_based.py @@ -569,19 +569,30 @@ def test_slitpore(self): system.part.add(pos=[0., 0., 0.], type=0) x = self.box_l / 2.0 + d = 1 - np.sqrt(2) / 2 parameters = [ ([x, x, 1.], -4., [0., 0., -1.]), # outside channel ([x, x, 15.], 5., [-1., 0., 0.]), # inside channel ([x, x, 5.], 0., [0., 0., 0.]), # on channel bottom surface + ([x - 5., x, 15.], 0., [0., 0., 0.]), # on channel side surface ([x + 5., x, 15.], 0., [0., 0., 0.]), # on channel side surface - ([x, x, 15.], 5., [-1., 0., 0.]), # within mouth - ([x, x, 25.], 0., [0., 0., 0.]), # on wall surface + ([x - 5. + 2 * d, x, 5. + 2 * d], 0., [0., 0., 0.]), # lower circle + ([x + 5. - 2 * d, x, 5. + 2 * d], 0., [0., 0., 0.]), # lower circle + ([x - 5. - 3 * d, x, 20. - 3 * d], 0., [0., 0., 0.]), # upper circle + ([x + 5. + 3 * d, x, 20. - 3 * d], 0., [0., 0., 0.]), # upper circle + ([1., x, 20.], 0., [0., 0., 0.]), # on inner wall surface + ([x, x, 25.], 0., [0., 0., 0.]), # on outer wall surface ([x, x, 27.], -2., [0., 0., 1.]), # outside wall ] for pos, ref_mindist, ref_force in parameters: system.part[0].pos = pos system.integrator.run(recalc_forces=True, steps=0) - self.assertEqual(slitpore_constraint.min_dist(), ref_mindist) + obs_mindist = slitpore_constraint.min_dist() + self.assertAlmostEqual(obs_mindist, ref_mindist, places=10) + if (ref_mindist == 0. and obs_mindist != 0.): + # force direction on a circle is not well-defined due to + # numerical instability + continue np.testing.assert_almost_equal( np.copy(slitpore_constraint.total_force()), ref_force, 10) @@ -792,7 +803,7 @@ def test_torus(self): # test summed forces on torus wall self.assertAlmostEqual( - -1.0 * torus_wall.total_force()[1], + torus_wall.total_force()[1], tests_common.lj_force( espressomd, cutoff=2.0, @@ -800,7 +811,7 @@ def test_torus(self): eps=1.0, sig=1.0, r=torus_constraint.min_dist()), - places=10) # minus for Newton's third law + places=10) # check whether total_summed_outer_normal_force is correct y_part2 = self.box_l / 2.0 + 2.0 * radius - part_offset diff --git a/testsuite/scripts/samples/CMakeLists.txt b/testsuite/scripts/samples/CMakeLists.txt index 2320b351b26..75178d45fa7 100644 --- a/testsuite/scripts/samples/CMakeLists.txt +++ b/testsuite/scripts/samples/CMakeLists.txt @@ -64,7 +64,7 @@ sample_test(FILE test_visualization_cellsystem.py) sample_test(FILE test_visualization_charged.py) foreach( var_constraint_shape - wall;sphere;ellipsoid;cylinder;spherocylinder;stomatocyte;simplepore;slitpore;hollowconicalfrustum + wall;sphere;ellipsoid;cylinder;spherocylinder;stomatocyte;simplepore;slitpore;torus;hollowconicalfrustum ) sample_test(FILE test_visualization_constraints.py SUFFIX ${var_constraint_shape})