Loading [MathJax]/extensions/TeX/AMSmath.js
ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages Concepts
Cylinder.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22#include <shapes/Cylinder.hpp>
23
24#include <utils/Vector.hpp>
25
26#include <cassert>
27#include <cmath>
28#include <utility>
29
30namespace Shapes {
31
32std::pair<double, double> Cylinder::dist_half_pore(double r, double z) const {
33 assert(z >= 0.0);
34 assert(r >= 0.0);
35
36 if (z >= m_half_length || r >= m_rad) {
37 /* Outside */
38 if (!m_open && z >= m_half_length && r < m_rad) {
39 /* Closest feature: cap */
40 return {0, -(z - m_half_length)};
41 }
42 if (z >= m_half_length && (m_open || r >= m_rad)) {
43 /* Closest feature: ring */
44 return {-(r - m_rad), -(z - m_half_length)};
45 }
46 /* Closest feature: cylinder */
47 return {-(r - m_rad), 0};
48 }
49 /* Inside */
50 if (!m_open && z >= m_half_length - m_rad &&
51 r < (z - (m_half_length - m_rad))) {
52 /* Closest feature: cap */
53 return {0, m_half_length - z};
54 }
55 /* Closest feature: cylinder */
56 return {m_rad - r, 0};
57}
58
59void Cylinder::calculate_dist(const Utils::Vector3d &pos, double &dist,
60 Utils::Vector3d &vec) const {
61 /* Coordinate transform to cylinder coords
62 with origin at m_center. */
63 Utils::Vector3d const c_dist = pos - m_center;
64 auto const z = e_z * c_dist;
65 auto const r_vec = c_dist - z * e_z;
66 auto const r = r_vec.norm();
67
68 /* If exactly on the axis, chose e_r orthogonal
69 to e_z. */
70 auto const e_r = (r == 0) ? e_r_axis : r_vec / r;
71
72 /* The pore has mirror symmetry in z with regard to
73 the center in the {r,z} system. We calculate always
74 for the z > 0 case, and flip the result if appropriate. */
75 auto [dr, dz] = dist_half_pore(r, std::abs(z));
76
77 double side = -1;
78 if (std::abs(z) >= m_half_length || r >= m_rad) /* outside */
79 side = 1;
80
81 if (z <= 0.0) {
82 dz *= -1;
83 }
84
85 dist = std::sqrt(dr * dr + dz * dz) * m_direction * side;
86 vec = -dr * e_r - dz * e_z;
87}
88} // namespace Shapes
Vector implementation and trait types for boost qvm interoperability.
bool m_open
whether to ignore bottom and top cap of cylinder
double m_half_length
Center of smoothing circle.
Utils::Vector3d e_r_axis
Alternative e_r for corner case.
Utils::Vector3d e_z
Unit vector in z direction.
void calculate_dist(const Utils::Vector3d &pos, double &dist, Utils::Vector3d &vec) const override
Definition Cylinder.cpp:59
double m_direction
direction -1: inside, +1 outside
Utils::Vector3d m_center
center of the cylinder.
std::pair< double, double > dist_half_pore(double r, double z) const
Definition Cylinder.cpp:32
T norm() const
Definition Vector.hpp:139