ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
Rhomboid.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/Rhomboid.hpp>
23
24#include <utils/Vector.hpp>
25
26#include <functional>
27
28namespace Shapes {
29void Rhomboid::calculate_dist(const Utils::Vector3d &pos, double &dist,
30 Utils::Vector3d &vec) const {
31
32 auto const le = std::less_equal<double>();
33 auto const ge = std::greater_equal<double>();
34 auto const lt = std::less<double>();
35 auto const gt = std::greater<double>();
36
37 // calculate vectors and scalars that are going to be used frequently
38
39 auto const axb = vector_product(m_a, m_b);
40 auto const bxc = vector_product(m_b, m_c);
41 auto const axc = vector_product(m_a, m_c);
42
43 auto const a_dot_bxc = m_a * bxc;
44 auto const b_dot_axc = m_b * axc;
45 auto const c_dot_axb = m_c * axb;
46
47 auto const dpos = pos - m_pos;
48
49 // compute distance from the rhomboid corners, edges and faces using linear
50 // combinations of the rhomboid edge vectors
51
52 auto const corner = [this, &vec, &dist, a = bxc / a_dot_bxc,
53 b = axc / b_dot_axc,
54 c = axb / c_dot_axb](auto op1, auto op2, auto op3,
55 Utils::Vector3d const &d) {
56 /* coefficients A, B, C tell whether ppos lies within a cone defined
57 * by pos and the adjacent edges */
58 auto const A = a * d;
59 auto const B = b * d;
60 auto const C = c * d;
61 if (op1(A, 0) and op2(B, 0) and op3(C, 0)) {
62 vec = d;
63 dist = m_direction * vec.norm();
64 return true;
65 }
66 return false;
67 };
68
69 if ( // check for cone at m_pos
70 corner(le, le, le, dpos) ||
71 // check for cone at m_pos+a
72 corner(ge, le, le, dpos - m_a) ||
73 // check for cone at m_pos+b
74 corner(le, ge, le, dpos - m_b) ||
75 // check for cone at m_pos+c
76 corner(le, le, ge, dpos - m_c) ||
77 // check for cone at m_pos+a+b
78 corner(ge, ge, le, dpos - m_a - m_b) ||
79 // check for cone at m_pos+a+c
80 corner(ge, le, ge, dpos - m_a - m_c) ||
81 // check for cone at m_pos+b+c
82 corner(le, ge, ge, dpos - m_b - m_c) ||
83 // check for cone at m_pos+a+b+c
84 corner(ge, ge, ge, dpos - m_a - m_b - m_c))
85 return;
86
87 auto const edge = [this, &vec, &dist](auto op1, auto op2,
88 Utils::Vector3d const &d,
89 Utils::Vector3d const &axis1,
90 double const dir1_dot_axis1,
91 Utils::Vector3d const &axis2,
92 double const dir2_dot_axis2,
93 Utils::Vector3d const &edge) {
94 auto const A = (d * axis1) / dir1_dot_axis1;
95 auto const B = (d * axis2) / dir2_dot_axis2;
96 if (op1(A, 0) and op2(B, 0)) {
97 auto const tmp = (d * edge) / edge.norm2();
98 vec = d - edge * tmp;
99 dist = m_direction * vec.norm();
100 return true;
101 }
102 return false;
103 };
104
105 if ( // check for prism at edge m_pos, a
106 edge(le, le, dpos, axc, b_dot_axc, axb, c_dot_axb, m_a) ||
107 // check for prism at edge m_pos, b
108 edge(le, le, dpos, bxc, a_dot_bxc, axb, c_dot_axb, m_b) ||
109 // check for prism at edge m_pos, c
110 edge(le, le, dpos, bxc, a_dot_bxc, axc, b_dot_axc, m_c) ||
111 // check for prism at edge m_pos+a, b
112 edge(ge, le, dpos - m_a, bxc, a_dot_bxc, axb, c_dot_axb, m_b) ||
113 // check for prism at edge m_pos+a, c
114 edge(ge, le, dpos - m_a, bxc, a_dot_bxc, axc, b_dot_axc, m_c) ||
115 // check for prism at edge m_pos+b+c, c
116 edge(le, ge, dpos - m_b - m_c, bxc, a_dot_bxc, axc, b_dot_axc, m_c) ||
117 // check for prism at edge m_pos+b+c, b
118 edge(le, ge, dpos - m_b - m_c, bxc, a_dot_bxc, axb, c_dot_axb, m_b) ||
119 // check for prism at edge m_pos+b+c, a
120 edge(ge, ge, dpos - m_b - m_c, axc, b_dot_axc, axb, c_dot_axb, m_a) ||
121 // check for prism at edge m_pos+a+b, a
122 edge(ge, le, dpos - m_a - m_b, axc, b_dot_axc, axb, c_dot_axb, m_a) ||
123 // check for prism at edge m_pos+a+b, c
124 edge(ge, ge, dpos - m_a - m_b, bxc, a_dot_bxc, axc, b_dot_axc, m_c) ||
125 // check for prism at edge m_pos+a+c, a
126 edge(le, ge, dpos - m_a - m_c, axc, b_dot_axc, axb, c_dot_axb, m_a) ||
127 // check for prism at edge m_pos+a+c, b
128 edge(ge, ge, dpos - m_a - m_c, bxc, a_dot_bxc, axb, c_dot_axb, m_b))
129 return;
130
131 auto const face_outside =
132 [this, &vec, &dist](auto op1, auto op2, Utils::Vector3d const &distance,
133 Utils::Vector3d const &axis,
134 double const dir_dot_axis, int sign) {
135 auto d = distance * axis;
136 if (op1(dir_dot_axis, 0)) {
137 d *= -1;
138 }
139 if (d >= 0) {
140 auto const tmp = axis.norm();
141 d /= tmp;
142 dist = d * m_direction;
143 if (op2(dir_dot_axis, 0)) {
144 d *= -1;
145 }
146 vec = (sign * d / tmp) * axis;
147 return true;
148 }
149 return false;
150 };
151
152 if ( // check for face with normal -axb
153 face_outside(gt, lt, dpos, axb, c_dot_axb, -1) ||
154 // calculate distance to face with normal axc
155 face_outside(gt, gt, dpos, axc, b_dot_axc, +1) ||
156 // calculate distance to face with normal -bxc
157 face_outside(gt, lt, dpos, bxc, a_dot_bxc, -1) ||
158 // calculate distance to face with normal axb
159 face_outside(lt, lt, dpos - m_a - m_b - m_c, axb, c_dot_axb, +1) ||
160 // calculate distance to face with normal -axc
161 face_outside(lt, gt, dpos - m_a - m_b - m_c, axc, b_dot_axc, -1) ||
162 // calculate distance to face with normal bxc
163 face_outside(lt, lt, dpos - m_a - m_b - m_c, bxc, a_dot_bxc, +1))
164 return;
165
166 // ppos lies within rhomboid.
167 // Find nearest wall for interaction (test all 6 possibilities).
168
169 // calculate distance to face with normal -axb
170 {
171 auto d = dpos * axb;
172 if (c_dot_axb > 0.0) {
173 d *= -1;
174 }
175 auto const tmp = axb.norm();
176 d /= tmp;
177 dist = d * m_direction;
178 if (c_dot_axb < 0.0) {
179 d *= -1;
180 }
181 vec = (-d / tmp) * axb;
182 }
183
184 auto const face_inside =
185 [this, &vec, &dist](auto op1, auto op2, Utils::Vector3d const &distance,
186 Utils::Vector3d const &axis,
187 double const dir_dot_axis, int sign) {
188 auto d = distance * axis;
189 if (op1(dir_dot_axis, 0)) {
190 d *= -1;
191 }
192 auto const tmp = axis.norm();
193 d /= tmp;
194 if (std::abs(d) < std::abs(dist)) {
195 dist = d * m_direction;
196 if (op2(dir_dot_axis, 0)) {
197 d *= -1;
198 }
199 vec = (sign * d / tmp) * axis;
200 }
201 };
202
203 // calculate distance to face with normal axc
204 face_inside(gt, gt, dpos, axc, b_dot_axc, +1);
205 // calculate distance to face with normal -bxc
206 face_inside(gt, lt, dpos, bxc, a_dot_bxc, -1);
207 // calculate distance to face with normal axb
208 face_inside(lt, lt, dpos - m_a - m_b - m_c, axb, c_dot_axb, +1);
209 // calculate distance to face with normal -axc
210 face_inside(lt, gt, dpos - m_a - m_b - m_c, axc, b_dot_axc, -1);
211 // calculate distance to face with normal bxc
212 face_inside(lt, lt, dpos - m_a - m_b - m_c, bxc, a_dot_bxc, +1);
213}
214
215} // namespace Shapes
Vector implementation and trait types for boost qvm interoperability.
void calculate_dist(const Utils::Vector3d &pos, double &dist, Utils::Vector3d &vec) const override
Definition Rhomboid.cpp:29
T norm() const
Definition Vector.hpp:138
__device__ void vector_product(float const *a, float const *b, float *out)