OpenLB 1.7
Loading...
Searching...
No Matches
boundarySimpleReflection3D.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2014 Robin Trunk, Mathias J. Krause
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program 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
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22 */
23
24#ifndef BOUNDARYSIMPLEREFLECTION3D_HH
25#define BOUNDARYSIMPLEREFLECTION3D_HH
26
27#include <set>
28#include "boundarySimpleReflection3D.h"
29
30namespace olb {
31
32template<typename T, template<typename U> class PARTICLETYPE>
33SimpleReflectBoundary3D<T, PARTICLETYPE>::SimpleReflectBoundary3D(T dT, SuperGeometry<T,3>& sg, std::set<int> materials) :
34 Boundary3D<T, PARTICLETYPE>(), _dT(dT), _sg(sg), _materials(materials)
35{
36 _matIter = _materials.begin();
37}
38
39template<typename T, template<typename U> class PARTICLETYPE>
40void SimpleReflectBoundary3D<T, PARTICLETYPE>::applyBoundary(typename std::deque<PARTICLETYPE<T> >::iterator& p, ParticleSystem3D<T, PARTICLETYPE>& psSys)
41{
42 std::vector<T> oldPos(3,T());
43 oldPos[0] = p->getPos()[0] - _dT * p->getVel()[0];
44 oldPos[1] = p->getPos()[1] - _dT * p->getVel()[1];
45 oldPos[2] = p->getPos()[2] - _dT * p->getVel()[2];
46
47 std::vector<T> line(3,T());
48 line[0] = p->getPos()[0] - oldPos[0];
49 line[1] = p->getPos()[1] - oldPos[1];
50 line[2] = p->getPos()[2] - oldPos[2];
51
52 std::vector<T> reflection(3, T());
53 int latticeR[4] = { 0,0,0,0 };
54 bool outer = false;
55
56 // get material number (Trigger)
57 // TODO: take particle radius into account for collision detection
58 latticeR[0] = p->getCuboid();
59 _sg.getCuboidGeometry().get(latticeR[0]).getLatticeR(&(latticeR[1]),&p->getPos()[0]);
60 int mat = _sg.get(latticeR);
61
62 for (_matIter = _materials.begin(); _matIter != _materials.end(); _matIter++) {
63 if (mat == *_matIter) {
64 // compute discrete normal
65 // TODO: rearrange the if to make computation more efficient
66 std::vector<int> normal(3, T());
67 if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2], latticeR[3]) == 1) {
68 normal[0] = -1;
69 }
70 if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2], latticeR[3]) == 1) {
71 normal[0] = 1;
72 }
73 if (_sg.get(latticeR[0], latticeR[1], latticeR[2] - 1, latticeR[3]) == 1) {
74 normal[1] = -1;
75 }
76 if (_sg.get(latticeR[0], latticeR[1], latticeR[2] + 1, latticeR[3]) == 1) {
77 normal[1] = 1;
78 }
79 if (_sg.get(latticeR[0], latticeR[1], latticeR[2], latticeR[3] - 1) == 1) {
80 normal[2] = -1;
81 }
82 if (_sg.get(latticeR[0], latticeR[1], latticeR[2], latticeR[3] + 1) == 1) {
83 normal[2] = 1;
84 }
85
86 if (normal[0]==0 && normal[1]==0 && normal[2]==0) {
87 outer = true;
88 // check for outer edge
89 if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] - 1, latticeR[3]) == 1) {
90 normal[0] = -1;
91 normal[1] = -1;
92 }
93 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] + 1, latticeR[3]) == 1) {
94 normal[0] = -1;
95 normal[1] = 1;
96 }
97 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] - 1, latticeR[3]) == 1) {
98 normal[0] = 1;
99 normal[1] = -1;
100 }
101 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] + 1, latticeR[3]) == 1) {
102 normal[0] = 1;
103 normal[1] = 1;
104 }
105 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2], latticeR[3] - 1) == 1) {
106 normal[0] = -1;
107 normal[2] = -1;
108 }
109 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2], latticeR[3] + 1) == 1) {
110 normal[0] = -1;
111 normal[2] = 1;
112 }
113 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2], latticeR[3] - 1) == 1) {
114 normal[0] = 1;
115 normal[2] = -1;
116 }
117 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2], latticeR[3] + 1) == 1) {
118 normal[0] = 1;
119 normal[2] = 1;
120 }
121 else if (_sg.get(latticeR[0], latticeR[1], latticeR[2] - 1, latticeR[3] - 1) == 1) {
122 normal[1] = -1;
123 normal[2] = -1;
124 }
125 else if (_sg.get(latticeR[0], latticeR[1], latticeR[2] - 1, latticeR[3] + 1) == 1) {
126 normal[1] = -1;
127 normal[2] = 1;
128 }
129 else if (_sg.get(latticeR[0], latticeR[1], latticeR[2] + 1, latticeR[3] - 1) == 1) {
130 normal[1] = +1;
131 normal[2] = -1;
132 }
133 else if (_sg.get(latticeR[0], latticeR[1], latticeR[2] + 1, latticeR[3] + 1) == 1) {
134 normal[1] = +1;
135 normal[2] = 1;
136 }
137 // check for outer corner
138 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] - 1, latticeR[3] - 1) == 1) {
139 normal[0] = -1;
140 normal[1] = -1;
141 normal[2] = -1;
142 }
143 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] - 1, latticeR[3] + 1) == 1) {
144 normal[0] = -1;
145 normal[1] = -1;
146 normal[2] = 1;
147 }
148 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] + 1, latticeR[3] - 1) == 1) {
149 normal[0] = -1;
150 normal[1] = +1;
151 normal[2] = -1;
152 }
153 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] + 1, latticeR[3] + 1) == 1) {
154 normal[0] = -1;
155 normal[1] = +1;
156 normal[2] = 1;
157 }
158 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] - 1, latticeR[3] - 1) == 1) {
159 normal[0] = +1;
160 normal[1] = -1;
161 normal[2] = -1;
162 }
163 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] - 1, latticeR[3] + 1) == 1) {
164 normal[0] = +1;
165 normal[1] = -1;
166 normal[2] = 1;
167 }
168 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] + 1, latticeR[3] - 1) == 1) {
169 normal[0] = +1;
170 normal[1] = +1;
171 normal[2] = -1;
172 }
173 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] + 1, latticeR[3] + 1) == 1) {
174 normal[0] = +1;
175 normal[1] = +1;
176 normal[2] = 1;
177 }
178 // Error since normal is still zero
179 else {
180 std::cout << "----->>>>> ERROR Normal is ZERO" << std::endl;
181 outer=false;
182 }
183 }
184
185 T prod = line[0]*normal[0] + line[1]*normal[1] + line[2]*normal[2];
186 // if particle has obtuse angle to normal do nothing
187 if (!outer) {
188 if (prod >= 0) {
189 return;
190 }
191 }
192
193 // if particle has acute angle to normal reverse normal component of velocity
194 T invNormNormal = 1. / util::norm2(normal);
195 reflection[0] = line[0] - 2.*prod*normal[0]*invNormNormal;
196 reflection[1] = line[1] - 2.*prod*normal[1]*invNormNormal;
197 reflection[2] = line[2] - 2.*prod*normal[2]*invNormNormal;
198
199 // compute new velocity vector
200 T vel = util::norm(p->getVel());
201 reflection = util::normalize(reflection);
202 reflection[0] = reflection[0]*vel;
203 reflection[1] = reflection[1]*vel;
204 reflection[2] = reflection[2]*vel;
205
206 // TODO: take distance to impact point into account, not just 0.5
207 oldPos[0] += 0.5 * (line[0] + reflection[0]) * _dT;
208 oldPos[1] += 0.5 * (line[1] + reflection[1]) * _dT;
209 oldPos[2] += 0.5 * (line[2] + reflection[2]) * _dT;
210
211 p->getVel() = reflection;
212 p->setPos(oldPos);
213
214 return;
215 }
216 }
217}
218
219
220template<typename T, template<typename U> class PARTICLETYPE>
221ReflectBoundary3D<T, PARTICLETYPE>::ReflectBoundary3D(T dT, SuperGeometry<T,3>& sg, std::set<int> materials) :
222 Boundary3D<T, PARTICLETYPE>(), _dT(dT), _sg(sg), _materials(materials)
223{
224 _matIter = _materials.begin();
225}
226
227template<typename T, template<typename U> class PARTICLETYPE>
228Vector<T,3> ReflectBoundary3D<T, PARTICLETYPE>::computeDiscreteNormal(int latticeR[])
229{
230 Vector<T,3> normal {0,0,0};
231 if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2], latticeR[3]) == 1) {
232 normal[0] = -1;
233 }
234 if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2], latticeR[3]) == 1) {
235 normal[0] = 1;
236 }
237 if (_sg.get(latticeR[0], latticeR[1], latticeR[2] - 1, latticeR[3]) == 1) {
238 normal[1] = -1;
239 }
240 if (_sg.get(latticeR[0], latticeR[1], latticeR[2] + 1, latticeR[3]) == 1) {
241 normal[1] = 1;
242 }
243 if (_sg.get(latticeR[0], latticeR[1], latticeR[2], latticeR[3] - 1) == 1) {
244 normal[2] = -1;
245 }
246 if (_sg.get(latticeR[0], latticeR[1], latticeR[2], latticeR[3] + 1) == 1) {
247 normal[2] = 1;
248 }
249
250 if (normal[0]==0 && normal[1]==0 && normal[2]==0) {
251 // check for outer edge. The normals to the planes forming the edge need to be pushed back to normal.
252 if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] - 1, latticeR[3]) == 1) {
253 normal[0] = -1;
254 normal[1] = -1;
255 }
256 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] + 1, latticeR[3]) == 1) {
257 normal[0] = -1;
258 normal[1] = 1;
259 }
260 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] - 1, latticeR[3]) == 1) {
261 normal[0] = 1;
262 normal[1] = -1;
263 }
264 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] + 1, latticeR[3]) == 1) {
265 normal[0] = 1;
266 normal[1] = 1;
267 }
268 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2], latticeR[3] - 1) == 1) {
269 normal[0] = -1;
270 normal[2] = -1;
271 }
272 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2], latticeR[3] + 1) == 1) {
273 normal[0] = -1;
274 normal[2] = 1;
275 }
276 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2], latticeR[3] - 1) == 1) {
277 normal[0] = 1;
278 normal[2] = -1;
279 }
280 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2], latticeR[3] + 1) == 1) {
281 normal[0] = 1;
282 normal[2] = 1;
283 }
284 else if (_sg.get(latticeR[0], latticeR[1], latticeR[2] - 1, latticeR[3] - 1) == 1) {
285 normal[1] = -1;
286 normal[2] = -1;
287 }
288 else if (_sg.get(latticeR[0], latticeR[1], latticeR[2] - 1, latticeR[3] + 1) == 1) {
289 normal[1] = -1;
290 normal[2] = 1;
291 }
292 else if (_sg.get(latticeR[0], latticeR[1], latticeR[2] + 1, latticeR[3] - 1) == 1) {
293 normal[1] = +1;
294 normal[2] = -1;
295 }
296 else if (_sg.get(latticeR[0], latticeR[1], latticeR[2] + 1, latticeR[3] + 1) == 1) {
297 normal[1] = +1;
298 normal[2] = 1;
299 }
300 // check for outer corner. The normals to the planes forming the corner need to be pushed back to normal.
301 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] - 1, latticeR[3] - 1) == 1) {
302 normal[0] = -1;
303 normal[1] = -1;
304 normal[2] = -1;
305 }
306 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] - 1, latticeR[3] + 1) == 1) {
307 normal[0] = -1;
308 normal[1] = -1;
309 normal[2] = 1;
310 }
311 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] + 1, latticeR[3] - 1) == 1) {
312 normal[0] = -1;
313 normal[1] = +1;
314 normal[2] = -1;
315 }
316 else if (_sg.get(latticeR[0], latticeR[1] - 1, latticeR[2] + 1, latticeR[3] + 1) == 1) {
317 normal[0] = -1;
318 normal[1] = +1;
319 normal[2] = 1;
320 }
321 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] - 1, latticeR[3] - 1) == 1) {
322 normal[0] = +1;
323 normal[1] = -1;
324 normal[2] = -1;
325 }
326 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] - 1, latticeR[3] + 1) == 1) {
327 normal[0] = +1;
328 normal[1] = -1;
329 normal[2] = 1;
330 }
331 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] + 1, latticeR[3] - 1) == 1) {
332 normal[0] = +1;
333 normal[1] = +1;
334 normal[2] = -1;
335 }
336 else if (_sg.get(latticeR[0], latticeR[1] + 1, latticeR[2] + 1, latticeR[3] + 1) == 1) {
337 normal[0] = +1;
338 normal[1] = +1;
339 normal[2] = 1;
340 }
341 // Error since normal[0] is still zero
342 else {
343 throw std::out_of_range("----->>>>> ERROR in ReflectionBoundary3D::computeDiscreteNormals(): Normal is ZERO");
344 }
345 }
346
347 return normal;
348}
349
350
351template<typename T, template<typename U> class PARTICLETYPE>
352void ReflectBoundary3D<T, PARTICLETYPE>::applyBoundary(typename std::deque<PARTICLETYPE<T> >::iterator& p, ParticleSystem3D<T, PARTICLETYPE>& psSys)
353{
354 //
355 int latticeR[4] = { 0,0,0,0 };
356 latticeR[0] = p->getCuboid();
357 _sg.getCuboidGeometry().get(latticeR[0]).getLatticeR(&(latticeR[1]),&p->getPos()[0]);
358
359 int mat = _sg.get(latticeR);
360 for (_matIter = _materials.begin(); _matIter != _materials.end(); _matIter++) {
361 if (mat == *_matIter) {
362 // Computing the discrete normal(s).
363 Vector<T,3> normal { computeDiscreteNormal(latticeR) };
364 Vector<T,3> P1 { p->getPos() };
365 Vector<T,3> R { _sg.getPhysR({latticeR[0], latticeR[1], latticeR[2], latticeR[3]}) };
366
367 std::vector<T> P1New { p->getPos() };
368 std::vector<T> vNew { p->getVel() };
369
370 for (int iD=0; iD<3; iD++) {
371 if (normal[iD]*(R[iD] - P1[iD]) > 0) {
372 P1New[iD] = 2.*R[iD] - P1[iD];
373 vNew[iD] = -vNew[iD];
374 }
375 }
376 p->setPos(P1New);
377 p->setVel(vNew);
378
379 return;
380 }
381 }
382}
383
384
385}
386#endif
olb::Boundary3D
Prototype for all particle boundaries.
Definition boundary3D.h:43
olb::ReflectBoundary3D::applyBoundary
void applyBoundary(typename std::deque< PARTICLETYPE< T > >::iterator &p, ParticleSystem3D< T, PARTICLETYPE > &psSys) override
Definition boundarySimpleReflection3D.hh:352
olb::ReflectBoundary3D::ReflectBoundary3D
ReflectBoundary3D(T dT, SuperGeometry< T, 3 > &sg, std::set< int > materials)
Definition boundarySimpleReflection3D.hh:221
olb::SimpleReflectBoundary3D::applyBoundary
void applyBoundary(typename std::deque< PARTICLETYPE< T > >::iterator &p, ParticleSystem3D< T, PARTICLETYPE > &psSys) override
Definition boundarySimpleReflection3D.hh:40
olb::SimpleReflectBoundary3D::SimpleReflectBoundary3D
SimpleReflectBoundary3D(T dT, SuperGeometry< T, 3 > &sg, std::set< int > materials)
Definition boundarySimpleReflection3D.hh:33
olb::SuperGeometry
Representation of a statistic for a parallel 2D geometry.
Definition superGeometry.h:70
olb::Vector
Plain old scalar vector.
Definition vector.h:47
olb::util::norm
T norm(const std::vector< T > &a)
l2 norm of a vector of arbitrary length
Definition vectorHelpers.h:146
olb::util::norm2
T norm2(const std::vector< T > &a)
l2 norm to the power of 2 of a vector of arbitrary length
Definition vectorHelpers.h:158
olb::util::normalize
Vector< T, D > normalize(const Vector< T, D > &a)
Definition vectorHelpers.h:193
olb
Top level namespace for all of OpenLB.
Definition boundaryPostProcessors2D.h:34
Generated on Tue Mar 5 2024 09:57:51 for OpenLB by doxygen 1.10.0