From 30ce98c9e0289b6394e2b9f5cf654bc7071c9f1b Mon Sep 17 00:00:00 2001 From: Michael Lahnert Date: Tue, 20 Nov 2018 12:09:03 +0100 Subject: [PATCH 1/3] Fix Simple Pore --- src/core/shapes/SimplePore.cpp | 129 ++++++++++++++++++++------------- 1 file changed, 80 insertions(+), 49 deletions(-) diff --git a/src/core/shapes/SimplePore.cpp b/src/core/shapes/SimplePore.cpp index 796810cf9c3..96b5f3ba0b9 100644 --- a/src/core/shapes/SimplePore.cpp +++ b/src/core/shapes/SimplePore.cpp @@ -30,56 +30,87 @@ namespace Shapes { * * @returns The distance vector from the surface in the cylinder system. */ -std::pair SimplePore::dist_half_pore(double r, double z) const { - assert(z >= 0.0); - assert(r >= 0.0); - - if (z <= c_z) { - /* Cylinder section */ - return {m_rad - r, 0}; - } else if ((r >= c_r) && (z >= m_half_length)) { - /* Wall section */ - return {0, m_half_length - z}; - } else { - /* Smoothing area */ - /* Vector to center of torus segment */ - auto const dr = c_r - r; - auto const dz = c_z - z; - - /* Rescale to surface */ - auto const d = std::sqrt(dr * dr + dz * dz); - auto const fac = (d - m_smoothing_rad) / d; - - return {fac * dr, fac * dz}; - } -} - -void SimplePore::calculate_dist(const Vector3d &pos, double *dist, - double *vec) const { - /* Coordinate transform to cylinder coords - with origin at m_center. */ - Vector3d const c_dist = pos - m_center; - auto const z = e_z * c_dist; - auto const r_vec = c_dist - z * e_z; - auto const r = r_vec.norm(); - - /* If exactly on the axis, chose e_r orthogonal - to e_z. */ - auto const e_r = (r == 0) ? e_r_axis : r_vec / r; - - /* The pore has mirror symmetry in z with regard to - the center in the {r,z} system. We calculate always - for the z > 0 case, and flip the result if appropriate. */ - double dr, dz; - std::tie(dr, dz) = dist_half_pore(r, std::abs(z)); - - if (z <= 0.0) { - dz *= -1; + std::pair SimplePore::dist_half_pore(double r, double z) const { + assert(z >= 0.0); + assert(r >= 0.0); + + /* + * We have to find the line that splits area 1 (r determines distance) from + * area 2 (z determines distance) inside pore. In area 3 we have to consider + * z and r to determine the distance. + * + * | x + * | 2 x + * | x + * _|_ x 1 + * \| ^ r + * 3 |-------| + * | z <- + */ + + if ((z <= c_z) && (r <= (c_z + c_r - z))) { + /* Cylinder section, inner */ + return {m_rad - r, 0}; + } + else if (((z >= c_z) && (r >= c_r)) || + ((z <= c_z) && (r > (c_z + c_r - z)))) { + /* Wall section and outer cylinder */ + return {0, m_half_length - z}; + } else { + /* Smoothing area */ + /* Vector to center of torus segment */ + auto const dr = c_r - r; + auto const dz = c_z - z; + + /* Rescale to surface */ + auto const d = std::sqrt(dr * dr + dz * dz); + auto const fac = (d - m_smoothing_rad) / d; + + return {fac * dr, fac * dz}; + } } - *dist = std::sqrt(dr * dr + dz * dz); - for (int i = 0; i < 3; i++) { - vec[i] = -dr * e_r[i] + -dz * e_z[i]; + void SimplePore::calculate_dist(const Vector3d &pos, double *dist, + double *vec) const { + /* Coordinate transform to cylinder coords + with origin at m_center. */ + Vector3d const c_dist = pos - m_center; + auto const z = e_z * c_dist; + auto const r_vec = c_dist - z * e_z; + auto const r = r_vec.norm(); + + /* If exactly on the axis, chose e_r orthogonal + to e_z. */ + auto const e_r = (r == 0) ? e_r_axis : r_vec / r; + + /* The pore has mirror symmetry in z with regard to + the center in the {r,z} system. We calculate always + for the z > 0 case, and flip the result if appropriate. */ + double dr, dz; + std::tie(dr, dz) = dist_half_pore(r, std::abs(z)); + + double side = -1; + if (((dz == 0) && (r <= m_rad)) || // cylinder section + ((dr == 0) && (std::abs(z) > m_half_length))) {// || + side = 1; + } else { + // smoothing area + if (std::abs(z) >= c_z) { + double angle = std::asin((std::abs(z) - c_z) / m_smoothing_rad); + double dist_offset = m_smoothing_rad - (std::cos(angle) * m_smoothing_rad); + if (m_half_length < std::abs(z) || r <= (m_rad + dist_offset)) { + side = 1; + } + } + } + + if (z <= 0.0) { + dz *= -1; + } + + *dist = std::sqrt(dr * dr + dz * dz) * side; + for (int i = 0; i < 3; i++) { + vec[i] = -dr * e_r[i] + -dz * e_z[i]; + } } -} } // namespace Shapes From 529e827b59b13670b21cab9147e562b058a8cd8a Mon Sep 17 00:00:00 2001 From: Michael Lahnert Date: Tue, 20 Nov 2018 12:23:17 +0100 Subject: [PATCH 2/3] Apply coding style patch --- src/core/shapes/SimplePore.cpp | 152 ++++++++++++++++----------------- 1 file changed, 76 insertions(+), 76 deletions(-) diff --git a/src/core/shapes/SimplePore.cpp b/src/core/shapes/SimplePore.cpp index 96b5f3ba0b9..3fc7f0bd80d 100644 --- a/src/core/shapes/SimplePore.cpp +++ b/src/core/shapes/SimplePore.cpp @@ -30,87 +30,87 @@ namespace Shapes { * * @returns The distance vector from the surface in the cylinder system. */ - std::pair SimplePore::dist_half_pore(double r, double z) const { - assert(z >= 0.0); - assert(r >= 0.0); - - /* - * We have to find the line that splits area 1 (r determines distance) from - * area 2 (z determines distance) inside pore. In area 3 we have to consider - * z and r to determine the distance. - * - * | x - * | 2 x - * | x - * _|_ x 1 - * \| ^ r - * 3 |-------| - * | z <- - */ - - if ((z <= c_z) && (r <= (c_z + c_r - z))) { - /* Cylinder section, inner */ - return {m_rad - r, 0}; - } - else if (((z >= c_z) && (r >= c_r)) || +std::pair SimplePore::dist_half_pore(double r, double z) const { + assert(z >= 0.0); + assert(r >= 0.0); + + /* + * We have to find the line that splits area 1 (r determines distance) from + * area 2 (z determines distance) inside pore. In area 3 we have to consider + * z and r to determine the distance. + * + * | x + * | 2 x + * | x + * _|_ x 1 + * \| ^ r + * 3 |-------| + * | z <- + */ + + if ((z <= c_z) && (r <= (c_z + c_r - z))) { + /* Cylinder section, inner */ + return {m_rad - r, 0}; + } else if (((z >= c_z) && (r >= c_r)) || ((z <= c_z) && (r > (c_z + c_r - z)))) { - /* Wall section and outer cylinder */ - return {0, m_half_length - z}; - } else { - /* Smoothing area */ - /* Vector to center of torus segment */ - auto const dr = c_r - r; - auto const dz = c_z - z; - - /* Rescale to surface */ - auto const d = std::sqrt(dr * dr + dz * dz); - auto const fac = (d - m_smoothing_rad) / d; - - return {fac * dr, fac * dz}; - } + /* Wall section and outer cylinder */ + return {0, m_half_length - z}; + } else { + /* Smoothing area */ + /* Vector to center of torus segment */ + auto const dr = c_r - r; + auto const dz = c_z - z; + + /* Rescale to surface */ + auto const d = std::sqrt(dr * dr + dz * dz); + auto const fac = (d - m_smoothing_rad) / d; + + return {fac * dr, fac * dz}; } - - void SimplePore::calculate_dist(const Vector3d &pos, double *dist, - double *vec) const { - /* Coordinate transform to cylinder coords - with origin at m_center. */ - Vector3d const c_dist = pos - m_center; - auto const z = e_z * c_dist; - auto const r_vec = c_dist - z * e_z; - auto const r = r_vec.norm(); - - /* If exactly on the axis, chose e_r orthogonal - to e_z. */ - auto const e_r = (r == 0) ? e_r_axis : r_vec / r; - - /* The pore has mirror symmetry in z with regard to - the center in the {r,z} system. We calculate always - for the z > 0 case, and flip the result if appropriate. */ - double dr, dz; - std::tie(dr, dz) = dist_half_pore(r, std::abs(z)); - - double side = -1; - if (((dz == 0) && (r <= m_rad)) || // cylinder section - ((dr == 0) && (std::abs(z) > m_half_length))) {// || - side = 1; - } else { - // smoothing area - if (std::abs(z) >= c_z) { - double angle = std::asin((std::abs(z) - c_z) / m_smoothing_rad); - double dist_offset = m_smoothing_rad - (std::cos(angle) * m_smoothing_rad); - if (m_half_length < std::abs(z) || r <= (m_rad + dist_offset)) { - side = 1; - } +} + +void SimplePore::calculate_dist(const Vector3d &pos, double *dist, + double *vec) const { + /* Coordinate transform to cylinder coords + with origin at m_center. */ + Vector3d const c_dist = pos - m_center; + auto const z = e_z * c_dist; + auto const r_vec = c_dist - z * e_z; + auto const r = r_vec.norm(); + + /* If exactly on the axis, chose e_r orthogonal + to e_z. */ + auto const e_r = (r == 0) ? e_r_axis : r_vec / r; + + /* The pore has mirror symmetry in z with regard to + the center in the {r,z} system. We calculate always + for the z > 0 case, and flip the result if appropriate. */ + double dr, dz; + std::tie(dr, dz) = dist_half_pore(r, std::abs(z)); + + double side = -1; + if (((dz == 0) && (r <= m_rad)) || // cylinder section + ((dr == 0) && (std::abs(z) > m_half_length))) { // || + side = 1; + } else { + // smoothing area + if (std::abs(z) >= c_z) { + double angle = std::asin((std::abs(z) - c_z) / m_smoothing_rad); + double dist_offset = + m_smoothing_rad - (std::cos(angle) * m_smoothing_rad); + if (m_half_length < std::abs(z) || r <= (m_rad + dist_offset)) { + side = 1; } } + } - if (z <= 0.0) { - dz *= -1; - } + if (z <= 0.0) { + dz *= -1; + } - *dist = std::sqrt(dr * dr + dz * dz) * side; - for (int i = 0; i < 3; i++) { - vec[i] = -dr * e_r[i] + -dz * e_z[i]; - } + *dist = std::sqrt(dr * dr + dz * dz) * side; + for (int i = 0; i < 3; i++) { + vec[i] = -dr * e_r[i] + -dz * e_z[i]; } +} } // namespace Shapes From d48fc009e7f71e4c19ae1465a642b59e135f446d Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 30 Nov 2018 17:21:36 +0100 Subject: [PATCH 3/3] testsuite: Added test for orientation of SimplePore. --- testsuite/python/simple_pore.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/testsuite/python/simple_pore.py b/testsuite/python/simple_pore.py index d33c85e1696..18cf31c05c0 100644 --- a/testsuite/python/simple_pore.py +++ b/testsuite/python/simple_pore.py @@ -31,7 +31,17 @@ "Features not available, skipping test!") class SimplePoreConstraint(ut.TestCase): - def test(self): + def test_orientation(self): + pore = SimplePore( + axis=[1., 0., 0.], radius=2., smoothing_radius=.1, length=2., center=[5., 5., 5.]) + + d, _ = pore.call_method("calc_distance", position=[.0, .0, .0]) + self.assertGreater(d, 0.) + + d, _ = pore.call_method("calc_distance", position=[5., 5., .0]) + self.assertLess(d, 0.) + + def test_stability(self): s = espressomd.System(box_l=[1.0, 1.0, 1.0]) s.seed = s.cell_system.get_state()['n_nodes'] * [1234] box_yz = 15.