ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
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:138