OpenLB 1.7
Loading...
Searching...
No Matches
particleSystem3D.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2016 Thomas Henn, Davide Dapelo
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 PARTICLESYSTEM_3D_HH
25#define PARTICLESYSTEM_3D_HH
26
27#include <iostream>
28#include <cstring>
29#include <algorithm>
30#include <sstream>
31#include <string>
32#include <unordered_map>
33#include <memory>
34#include <stdexcept>
35
36#include "particle3D.h"
37#include "particleSystem3D.h"
38#include "functors/analytical/frameChangeF3D.h"
39 //#include "../functors/frameChangeF3D.h" //check
40 //#include "../utilities/vectorHelpers.h" //check
41
42#ifndef M_PI
43#define M_PI 3.14159265358979323846
44#endif
45
46// For agglomeration functions
47
48namespace olb {
49
50 template<typename T, template<typename U> class PARTICLETYPE>
51 ParticleSystem3D<T, PARTICLETYPE>::ParticleSystem3D( int iGeometry,
52 SuperGeometry<T,3>& superGeometry)
53 : clout(std::cout, "ParticleSystem3D"),
54 _iGeometry(iGeometry),
55 _superGeometry(superGeometry),
56 _contactDetection(new ContactDetection<T, PARTICLETYPE>(*this)),
57 _sim(this)
58 {
59 }
60
61 template<typename T, template<typename U> class PARTICLETYPE>
62 ParticleSystem3D<T, PARTICLETYPE>::ParticleSystem3D(
63 const ParticleSystem3D<T, PARTICLETYPE>& pS)
64 : clout(std::cout, "ParticleSystem3D"),
65 _iGeometry(pS._iGeometry),
66 _superGeometry(pS._superGeometry),
67 _contactDetection(new ContactDetection<T, PARTICLETYPE>(*this)),
68 _sim(this),
69 _physPos(pS._physPos),
70 _physExtend(pS._physExtend)
71 {
72 _particles = pS._particles;
73 _forces = pS._forces;
74 }
75
76 template<typename T, template<typename U> class PARTICLETYPE>
77 ParticleSystem3D<T, PARTICLETYPE>::ParticleSystem3D(
78 ParticleSystem3D<T, PARTICLETYPE> && pS)
79 : clout(std::cout, "ParticleSystem3D"),
80 _iGeometry(pS._iGeometry),
81 _superGeometry(pS._superGeometry),
82 _contactDetection(new ContactDetection<T, PARTICLETYPE>(*this)),
83 _sim(this),
84 _physPos(pS._physPos),
85 _physExtend(pS._physExtend)
86 {
87 _particles = std::move(pS._particles);
88 _forces = std::move(pS._forces);
89 }
90
91 template<typename T, template<typename U> class PARTICLETYPE>
92 ContactDetection<T, PARTICLETYPE>* ParticleSystem3D<T, PARTICLETYPE>::getContactDetection()
93 {
94 return _contactDetection;
95 }
96
97 template<typename T, template<typename U> class PARTICLETYPE>
98 void ParticleSystem3D<T, PARTICLETYPE>::printDeep(std::string message)
99 {
100 std::deque<PARTICLETYPE<T>*> allParticles = this->getAllParticlesPointer();
101 for (auto p : allParticles) {
102 p->printDeep("");
103 }
104 }
105
106 template<typename T, template<typename U> class PARTICLETYPE>
107 void ParticleSystem3D<T, PARTICLETYPE>::setContactDetection(ContactDetection<T, PARTICLETYPE>& contactDetection)
108 {
109 delete _contactDetection;
110 _contactDetection = contactDetection.generate(*this);
111 }
112
113 template<typename T, template<typename U> class PARTICLETYPE>
114 void ParticleSystem3D<T, PARTICLETYPE>::setPosExt(std::vector<T> physPos,
115 std::vector<T> physExtend)
116 {
117 _physPos = physPos;
118 _physExtend = physExtend;
119 }
120
121 template<typename T, template<typename U> class PARTICLETYPE>
122 const std::vector<T>& ParticleSystem3D<T, PARTICLETYPE>::getPhysPos()
123 {
124 return _physPos;
125 }
126
127 template<typename T, template<typename U> class PARTICLETYPE>
128 const std::vector<T>& ParticleSystem3D<T, PARTICLETYPE>::getPhysExtend()
129 {
130 return _physExtend;
131 }
132
133 template<typename T, template<typename U> class PARTICLETYPE>
134 PARTICLETYPE<T>& ParticleSystem3D<T, PARTICLETYPE>::operator[](const int i)
135 {
136 if (i < (int) _particles.size()) {
137 return _particles[i];
138 }
139 else if (i < (int) (_particles.size() + _shadowParticles.size())) {
140 return _shadowParticles[i - _particles.size()];
141 }
142 else {
143 std::ostringstream os;
144 os << i;
145 std::string err = "Cannot access element: " + os.str()
146 + " to large";
147 throw std::runtime_error(err);
148 }
149 }
150
151 template<typename T, template<typename U> class PARTICLETYPE>
152 const PARTICLETYPE<T>& ParticleSystem3D<T, PARTICLETYPE>::operator[](
153 const int i) const
154 {
155 if ((unsigned)i < _particles.size()) {
156 return _particles[i];
157 }
158 else if (i - _particles.size() < _shadowParticles.size()) {
159 return _shadowParticles[i - _particles.size()];
160 }
161 else {
162 std::ostringstream os;
163 os << i;
164 std::string err = "Cannot access element: " + os.str()
165 + " to large";
166 throw std::runtime_error(err);
167 }
168 }
169
170 template<typename T, template<typename U> class PARTICLETYPE>
171 int ParticleSystem3D<T, PARTICLETYPE>::size()
172 {
173 return _particles.size();
174 }
175
176 template<typename T, template<typename U> class PARTICLETYPE>
177 int ParticleSystem3D<T, PARTICLETYPE>::sizeInclShadow() const
178 {
179 return _particles.size() + _shadowParticles.size();
180 }
181
182 template<typename T, template<typename U> class PARTICLETYPE>
183 int ParticleSystem3D<T, PARTICLETYPE>::numOfActiveParticles()
184 {
185 int activeParticles = 0;
186 for (auto p : _particles) {
187 if (p.getActive()) {
188 activeParticles++;
189 }
190 }
191 return activeParticles;
192 }
193 /*
194 template<typename T, template<typename U> class PARTICLETYPE>
195 std::map<T, int> ParticleSystem3D<T, PARTICLETYPE>::radiusDistribution()
196 {
197 std::map<T, int> radDistr;
198 typename std::map<T, int>::iterator it;
199 T rad = 0;
200 for (auto p : _particles) {
201 if (!p.getAggl()) {
202 rad = p.getRad();
203 if (radDistr.count(rad) > 0) {
204 it = radDistr.find(rad);
205 it->second = it->second + 1;
206 } else {
207 radDistr.insert(std::pair<T, int>(rad, 1));
208 }
209 }
210 }
211 return radDistr;
212 }
213 */
214 template<typename T, template<typename U> class PARTICLETYPE>
215 int ParticleSystem3D<T, PARTICLETYPE>::numOfForces()
216 {
217 return _forces.size();
218 }
219
220 template<typename T, template<typename U> class PARTICLETYPE>
221 void ParticleSystem3D<T, PARTICLETYPE>::addParticle(PARTICLETYPE<T> &p)
222 {
223 _particles.push_back(p);
224 }
225
226 template<typename T, template<typename U> class PARTICLETYPE>
227 void ParticleSystem3D<T, PARTICLETYPE>::clearParticles()
228 {
229 _particles.clear();
230 }
231
232 template<typename T, template<typename U> class PARTICLETYPE>
233 void ParticleSystem3D<T, PARTICLETYPE>::addShadowParticle(PARTICLETYPE<T> &p)
234 {
235 _shadowParticles.push_back(p);
236 }
237
238 template<typename T, template<typename U> class PARTICLETYPE>
239 void ParticleSystem3D<T, PARTICLETYPE>::addForce(
240 std::shared_ptr<Force3D<T, PARTICLETYPE> > pF)
241 {
242 _forces.push_back(pF);
243 }
244
245 template<typename T, template<typename U> class PARTICLETYPE>
246 void ParticleSystem3D<T, PARTICLETYPE>::addBoundary(
247 std::shared_ptr<Boundary3D<T, PARTICLETYPE> > pB)
248 {
249 _boundaries.push_back(pB);
250 }
251
252 template<typename T, template<typename U> class PARTICLETYPE>
253 void ParticleSystem3D<T, PARTICLETYPE>::addParticleOperation(
254 std::shared_ptr<ParticleOperation3D<T, PARTICLETYPE> > pO)
255 {
256 _particleOperations.push_back(pO);
257 }
258
259 template<typename T, template<typename U> class PARTICLETYPE>
260 template<typename DESCRIPTOR>
261 void ParticleSystem3D<T, PARTICLETYPE>::setVelToFluidVel(
262 SuperLatticeInterpPhysVelocity3D<T, DESCRIPTOR> & fVel)
263 {
264 for (auto& p : _particles) {
265 if (p.getActive()) {
266 fVel(&p.getVel()[0], &p.getPos()[0], p.getCuboid());
267 }
268 }
269 }
270
271 template<typename T, template<typename U> class PARTICLETYPE>
272 void ParticleSystem3D<T, PARTICLETYPE>::setVelToAnalyticalVel(
273 AnalyticalConst3D<T, T>& aVel)
274 {
275 for (auto& p : _particles) {
276 if (p.getActive()) {
277 std::vector<T> vel(3, T());
278 aVel(&p.getVel()[0], &p.getPos()[0]);
279 }
280 }
281 }
282
283 template<typename T, template<typename U> class PARTICLETYPE>
284 void ParticleSystem3D<T, PARTICLETYPE>::computeForce()
285 {
286 typename std::deque<PARTICLETYPE<T> >::iterator p;
287 int pInt = 0;
288 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
289 if (p->getActive()) {
290 p->resetForce();
291 for (auto f : _forces) {
292 f->applyForce(p, pInt, *this);
293 }
294 }
295 }
296 }
297
298 template<typename T, template<typename U> class PARTICLETYPE>
299 void ParticleSystem3D<T, PARTICLETYPE>::eraseInactiveParticles()
300 {
301 std::deque<PARTICLETYPE<T> > activeParticles;
302 for (auto p = _particles.begin(); p != _particles.end(); ++p) {
303 if (p-> getActive()) {
304 PARTICLETYPE<T> particle = *p;
305 activeParticles.push_back(particle);
306 }
307 }
308 _particles = std::move(activeParticles);
309 }
310
311 template<typename T, template<typename U> class PARTICLETYPE>
312 void ParticleSystem3D<T, PARTICLETYPE>::computeBoundary()
313 {
314 typename std::deque<PARTICLETYPE<T> >::iterator p;
315 for (auto f : _boundaries) {
316 for (p = _particles.begin(); p != _particles.end(); ++p) {
317 if (p->getActive()) {
318 f->applyBoundary(p, *this);
319 }
320 }
321 }
322 }
323
324 template<typename T, template<typename U> class PARTICLETYPE>
325 void ParticleSystem3D<T, PARTICLETYPE>::computeParticleOperation()
326 {
327 typename std::deque<PARTICLETYPE<T> >::iterator p;
328 for (auto o : _particleOperations) {
329 for (p = _particles.begin(); p != _particles.end(); ++p) {
330 if (p->getActive()) {
331 o->applyParticleOperation(p, *this);
332 }
333 }
334 }
335 }
336
337 template<typename T, template<typename U> class PARTICLETYPE>
338 void ParticleSystem3D<T, PARTICLETYPE>::simulate(T dT, bool scale)
339 {
340 _sim.simulate(dT, scale);
341 }
342
343 template<typename T, template<typename U> class PARTICLETYPE>
344 void ParticleSystem3D<T, PARTICLETYPE>::simulateWithTwoWayCoupling_Mathias ( T dT,
345 ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling,
346 BackCouplingModel<T,PARTICLETYPE>& backCoupling,
347 int material, int subSteps, bool scale )
348 {
349 _sim.simulateWithTwoWayCoupling_Mathias(dT, forwardCoupling, backCoupling, material, subSteps, scale);
350 }
351
352 template<typename T, template<typename U> class PARTICLETYPE>
353 void ParticleSystem3D<T, PARTICLETYPE>::simulateWithTwoWayCoupling_Davide ( T dT,
354 ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling,
355 BackCouplingModel<T,PARTICLETYPE>& backCoupling,
356 int material, int subSteps, bool scale )
357 {
358 _sim.simulateWithTwoWayCoupling_Davide(dT, forwardCoupling, backCoupling, material, subSteps, scale);
359 }
360
361 // multiple collision models
362 template<typename T, template<typename U> class MagneticParticle3D>
363 void ParticleSystem3D<T, MagneticParticle3D>::simulate(T dT, std::set<int> sActivity, bool scale)
364 {
365 _sim.simulate(dT, sActivity, scale);
366 }
367
368 template<typename T, template<typename U> class PARTICLETYPE>
369 int ParticleSystem3D<T, PARTICLETYPE>::countMaterial(int mat)
370 {
371 int num = 0;
372 int locLatCoords[3] = {0};
373 for (auto& p : _particles) {
374 _superGeometry.getCuboidGeometry().get(p.getCuboid()).getFloorLatticeR(
375 locLatCoords, &p.getPos()[0]);
376 const BlockGeometry<T,3>& bg = _superGeometry.getBlockGeometry(_superGeometry.getLoadBalancer().loc(p.getCuboid()));
377 int iX = locLatCoords[0];
378 int iY = locLatCoords[1];
379 int iZ = locLatCoords[2];
380 if (bg.get(iX, iY, iZ) == mat
381 || bg.get(iX, iY + 1, iZ) == mat
382 || bg.get(iX, iY, iZ + 1) == mat
383 || bg.get(iX, iY + 1, iZ + 1) == mat
384 || bg.get(iX + 1, iY, iZ) == mat
385 || bg.get(iX + 1, iY + 1, iZ) == mat
386 || bg.get(iX + 1, iY, iZ + 1) == mat
387 || bg.get(iX + 1, iY + 1, iZ + 1) == mat) {
388 num++;
389 }
390 }
391 return num;
392 }
393
394 template<typename T, template<typename U> class PARTICLETYPE>
395 bool ParticleSystem3D<T, PARTICLETYPE>::executeForwardCoupling(ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling)
396 {
397 for (auto& p : _particles) {
398 if (p.getActive()) {
399 if (! forwardCoupling(&p, _iGeometry) ) {
400 std::cout << "ERROR: ForwardCoupling functional failed on particle " << p.getID();
401 return false;
402 }
403 }
404 }
405 return true;
406 }
407
408 template<typename T, template<typename U> class PARTICLETYPE>
409 bool ParticleSystem3D<T, PARTICLETYPE>::executeBackwardCoupling(BackCouplingModel<T,PARTICLETYPE>& backCoupling, int material, int subSteps)
410 {
411 for (auto& p : _particles) {
412 if (p.getActive()) {
413 if (! backCoupling(&p, _iGeometry, material, subSteps) ) {
414 std::cout << "ERROR: BackwardCoupling functional failed on particle " << p.getID();
415 return false;
416 }
417 }
418 }
419 return true;
420 }
421
422 template<typename T, template<typename U> class PARTICLETYPE>
423 void ParticleSystem3D<T, PARTICLETYPE>::explicitEuler(T dT, bool scale)
424 {
425 //std::cout << "explicit Euler" << std::endl;
426
427 T maxDeltaR = _superGeometry.getCuboidGeometry().getMaxDeltaR();
428 T maxFactor = T();
429
430 for (auto& p : _particles) {
431
432 if (p.getActive()) {
433
434 for (int i = 0; i < 3; i++) {
435 // std::cout <<p.getInvEffectiveMass()<< std::endl;
436 p.getVel()[i] += p.getForce()[i] * p.getInvEffectiveMass() * dT;
437 //std::cout << p.getVel()[i] << std::endl;
438 p.getPos()[i] += p.getVel()[i] * dT;
439 // std::cout << p.getPos()[i] << std::endl;
440
441 // computation of direction depending maxFactor to scale velocity value
442 // to keep velocity small enough for simulation
443 if (util::fabs(p.getVel()[i]) > util::fabs(maxDeltaR / dT)) {
444 maxFactor = util::max(maxFactor, util::fabs(p.getVel()[i] / maxDeltaR * dT));
445 }
446 }
447 // scaling of velocity values
448 // if particles are too fast, e.g. material boundary can not work anymore
449 if (!util::nearZero(maxFactor) && scale) {
450 std::cout << "particle velocity is scaled because of reached limit"
451 << std::endl;
452 for (int i = 0; i < 3; i++) {
453 p.getPos()[i] -= p.getVel()[i] * dT; // position set back
454 p.getVel()[i] /= maxFactor; // scale velocity value
455 p.getPos()[i] += p.getVel()[i] * dT;
456 }
457 }
458
459
460 // if particles are too fast, e.g. material boundary can not work anymore
461 //#ifdef OLB_DEBUG
462 // if (p.getVel()[i] * dT > _superGeometry.getCuboidGeometry().getMaxDeltaR()) {
463 // std::cout << " PROBLEM: particle speed too high rel. to delta of "
464 // "lattice: "<< std::endl;
465 // std::cout << "p.getVel()[i]*dT: " << i <<" "<< p.getVel()[i] * dT;
466 // std::cout << "MaxDeltaR(): " <<
467 // _superGeometry.getCuboidGeometry().getMaxDeltaR() << std::endl;
468 // exit(-1);
469 // }
470 //#endif
471
472 }
473 }
474 }
475
476 template<>
477 void ParticleSystem3D<double, MagneticParticle3D>::explicitEuler(double dT, bool scale)
478 {
479 double maxDeltaR = _superGeometry.getCuboidGeometry().getMaxDeltaR();
480 double maxFactor = double();
481
482 for (auto& p : _particles) {
483
484 if (p.getActive()) {
485
486 if (p.getSActivity() == 3) {
487 continue;
488 }
489
490 for (int i = 0; i < 3; i++) {
491 p.getVel()[i] += p.getForce()[i] * p.getInvEffectiveMass() * dT;
492 p.getPos()[i] += p.getVel()[i] * dT;
493
494 // computation of direction depending maxFactor to scale velocity value
495 // to keep velocity small enough for simulation
496 if (util::fabs(p.getVel()[i]) > util::fabs(maxDeltaR / dT)) {
497 maxFactor = util::max(maxFactor, util::fabs(p.getVel()[i] / maxDeltaR * dT));
498 }
499 }
500 // scaling of velocity values
501 // if particles are too fast, e.g. material boundary can not work anymore
502 if ( !util::nearZero(maxFactor) && scale) {
503 std::cout << "particle velocity is scaled because of reached limit"
504 << std::endl;
505 for (int i = 0; i < 3; i++) {
506 p.getPos()[i] -= p.getVel()[i] * dT; // position set back
507 p.getVel()[i] /= maxFactor; // scale velocity value
508 p.getPos()[i] += p.getVel()[i] * dT;
509 }
510 }
511
512 // Change sActivity of particle in dependence of its position in the geometry
513 // if ((p.getPos()[0] > 0.00015) && (p.getSActivity() == 0)) {
514 // p.setSActivity(2);
515 // }
516
517 // }
518
519 // if particles are too fast, e.g. material boundary can not work anymore
520 //#ifdef OLB_DEBUG
521 // if (p.getVel()[i] * dT > _superGeometry.getCuboidGeometry().getMaxDeltaR()) {
522 // std::cout << " PROBLEM: particle speed too high rel. to delta of "
523 // "lattice: "<< std::endl;
524 // std::cout << "p.getVel()[i]*dT: " << i <<" "<< p.getVel()[i] * dT;
525 // std::cout << "MaxDeltaR(): " <<
526 // _superGeometry.getCuboidGeometry().getMaxDeltaR() << std::endl;
527 // exit(-1);
528 // }
529 //#endif
530
531 }
532 }
533 }
534
535 // multiple collision models
536 template<>
537 void ParticleSystem3D<double, MagneticParticle3D>::explicitEuler(double dT, std::set<int> sActivityOfParticle, bool scale)
538 {
539 double maxDeltaR = _superGeometry.getCuboidGeometry().getMaxDeltaR();
540 double maxFactor = double();
541
542 for (auto& p : _particles) {
543
544 if (p.getActive()) {
545
546 if (p.getSActivity() == 3) {
547 continue;
548 }
549
550 bool b = false;
551 for (auto sA : sActivityOfParticle) {
552 if (p.getSActivity() == sA) {
553 b = true;
554 break;
555 }
556 }
557 if (b == false) {
558 continue;
559 }
560
561 for (int i = 0; i < 3; i++) {
562 p.getVel()[i] += p.getForce()[i] * p.getInvEffectiveMass() * dT;
563 p.getPos()[i] += p.getVel()[i] * dT;
564
565 // computation of direction depending maxFactor to scale velocity value
566 // to keep velocity small enough for simulation
567 if (util::fabs(p.getVel()[i]) > util::fabs(maxDeltaR / dT)) {
568 maxFactor = util::max(maxFactor, util::fabs(p.getVel()[i] / maxDeltaR * dT));
569 }
570 }
571 // scaling of velocity values
572 // if particles are too fast, e.g. material boundary can not work anymore
573 if ( !util::nearZero(maxFactor) && scale) {
574 std::cout << "particle velocity is scaled because of reached limit"
575 << std::endl;
576 for (int i = 0; i < 3; i++) {
577 p.getPos()[i] -= p.getVel()[i] * dT; // position set back
578 p.getVel()[i] /= maxFactor; // scale velocity value
579 p.getPos()[i] += p.getVel()[i] * dT;
580 }
581 }
582
583 // Change sActivity of particle in dependence of its position in the geometry
584 // if ((p.getPos()[0] > 0.00015) && (p.getSActivity() == 0)) {
585 // p.setSActivity(2);
586 // }
587
588 // }
589
590 // if particles are too fast, e.g. material boundary can not work anymore
591 //#ifdef OLB_DEBUG
592 // if (p.getVel()[i] * dT > _superGeometry.getCuboidGeometry().getMaxDeltaR()) {
593 // std::cout << " PROBLEM: particle speed too high rel. to delta of "
594 // "lattice: "<< std::endl;
595 // std::cout << "p.getVel()[i]*dT: " << i <<" "<< p.getVel()[i] * dT;
596 // std::cout << "MaxDeltaR(): " <<
597 // _superGeometry.getCuboidGeometry().getMaxDeltaR() << std::endl;
598 // exit(-1);
599 // }
600 //#endif
601
602 }
603 }
604 }
605
606 //template<typename T, template<typename U> class PARTICLETYPE>
607 //void ParticleSystem3D<T, PARTICLETYPE>::integrateTorque(T dT)
608 //{
609 // for (auto& p : _particles) {
610 // if (p.getActive()) {
611 // for (int i = 0; i < 3; i++) {
612 // p.getAVel()[i] += p.getTorque()[i] * 1.
613 // / (2. / 5. * p.getMass() * util::pow(p.getRad(), 2)) * dT;
614 // }
615 // }
616 // }
617 //}
618
619 //template<typename T, template<typename U> class PARTICLETYPE>
620 //void ParticleSystem3D<T, PARTICLETYPE>::integrateTorqueMag(T dT) {
623 // for (auto& p : _particles) {
624 // if (p.getActive()) {
625 // Vector<T, 3> deltaAngle;
626 // T angle;
627 // T epsilon = std::numeric_limits<T>::epsilon();
628 // for (int i = 0; i < 3; i++) {
629 // // change orientation due to torque moments
630 // deltaAngle[i] = (5. * p.getTorque()[i] * dT * dT) / (2. * p.getMass() * util::pow(p.getRad(), 2));
631 // // apply change in angle to dMoment vector
632 // }
633 // angle = norm(deltaAngle);
634 // if (angle > epsilon) {
635 // //Vector<T, 3> axis(deltaAngle);
636 // // Vector<T, 3> axis(T(), T(), T(1));
637 // //axis.normalize();
638 // std::vector<T> null(3, T());
639 //
640 // //RotationRoundAxis3D<T, S> rotRAxis(p.getPos(), util::fromVector3(axis), angle);
641 // RotationRoundAxis3D<T, S> rotRAxis(null, util::fromVector3(deltaAngle), angle);
642 // T input[3] = {p.getMoment()[0], p.getMoment()[1], p.getMoment()[2]};
643 // Vector<T, 3> in(input);
644 // T output[3] = {T(), T(), T()};
645 // rotRAxis(output, input);
646 // Vector<T, 3> out(output);
654 // std::cout<< "|moment|_in: " << in.norm() << ", |moment|_out: " << out.norm()
655 // << ", | |in| - |out| |: " << util::fabs(in.norm() - out.norm())
656 // /*<< " Angle: " << angle[0] */<< std::endl;
657 // p.getMoment()[0] = output[0];
658 // p.getMoment()[1] = output[1];
659 // p.getMoment()[2] = output[2];
660 // }
661 // }
662 // }
663 //}
664
665 //template<typename T, template<typename U> class PARTICLETYPE>
666 //void ParticleSystem3D<T, PARTICLETYPE>::implicitEuler(T dT, AnalyticalF<3,T,T>& getvel) {
667 // _activeParticles = 0;
668 // for (auto& p : _particles) {
669 // if(p.getActive()) {
670 // std::vector<T> fVel = getvel(p._pos);
672 // std::vector<T> pos = p.getPos();
673 // T C = 6.* M_PI * p._rad * this->_converter.getCharNu() * this->_converter.getCharRho()* dT/p.getMass();
674 // for (int i = 0; i<3; i++) {
675 // p._vel[i] = (p._vel[i]+C*fVel[i]) / (1+C);
676 // p._pos[i] += p._vel[i] * dT;
677 // }
680 // checkActive(p);
682 // }
684 // }
685 //}
686
687 /*
688 template<typename T, template<typename U> class PARTICLETYPE>
689 void ParticleSystem3D<T, PARTICLETYPE>::predictorCorrector1(T dT)
690 {
691 std::vector<T> vel;
692 std::vector<T> pos;
693 std::vector<T> frc;
694 for (auto& p : _particles) {
695 if (p.getActive()) {
696 vel = p.getVel();
697 p.setVel(vel, 1);
698 pos = p.getPos();
699 p.setVel(pos, 2);
700 frc = p.getForce();
701 p.setForce(frc, 1);
702 for (int i = 0; i < 3; i++) {
703 vel[i] += p._force[i] / p.getMass() * dT;
704 pos[i] += vel[i] * dT;
705 }
706 p.setVel(vel);
707 p.setPos(pos);
708 }
709 }
710 }
711
712 template<typename T, template<typename U> class PARTICLETYPE>
713 void ParticleSystem3D<T, PARTICLETYPE>::predictorCorrector2(T dT)
714 {
715 std::vector<T> vel;
716 std::vector<T> pos;
717 std::vector<T> frc;
718 for (auto& p : _particles) {
719 if (p.getActive()) {
720 vel = p.getVel(1);
721 pos = p.getVel(2);
722 for (int i = 0; i < 3; i++) {
723 vel[i] += dT * .5 * (p.getForce()[i] + p.getForce(1)[i]) / p.getMass();
724 pos[i] += vel[i] * dT;
725 }
726 p.setVel(vel);
727 p.setPos(pos);
728 }
729 }
730 }
731 */
732
733 /*
734 template<typename T, template<typename U> class PARTICLETYPE>
735 void ParticleSystem3D<T, PARTICLETYPE>::adamBashforth4(T dT)
736 {
737 for (auto& p : _particles) {
738 if (p.getActive()) {
739 std::vector<T> vel = p.getVel();
740 std::vector<T> pos = p.getPos();
741 for (int i = 0; i < 3; i++) {
742 vel[i] += dT / p.getMas()
743 * (55. / 24. * p.getForce()[i] - 59. / 24. * p.getForce(1)[i]
744 + 37. / 24. * p.getForce(2)[i] - 9. / 24. * p.getForce(3)[i]);
745 }
746 p.rotAndSetVel(vel);
747 for (int i = 0; i < 3; i++) {
748 pos[i] += dT
749 * (55. / 24. * p.getVel()[i] - 59. / 24. * p.getVel(1)[i]
750 + 37. / 24. * p.getVel(2)[i] - 3. / 8. * p.getVel(3)[i]);
751 }
752 p.setPos(pos);
753 }
754 }
755 }
756 */
757
758 template<typename T, template<typename U> class PARTICLETYPE>
759 void ParticleSystem3D<T, PARTICLETYPE>::velocityVerlet1(T dT)
760 {
761 for (auto& p : _particles) {
762 if (p.getActive()) {
763 std::vector<T> pos = p.getPos();
764 for (int i = 0; i < 3; i++) {
765 pos[i] += p.getVel()[i] * dT
766 + .5 * p.getForce()[i] / p.getMass() * util::pow(dT, 2);
767 }
768 p.setPos(pos);
769 }
770 }
771 }
772
773 template<typename T, template<typename U> class PARTICLETYPE>
774 void ParticleSystem3D<T, PARTICLETYPE>::velocityVerlet2(T dT)
775 {
776 std::vector<T> frc;
777 std::vector<T> vel;
778 typename std::deque<PARTICLETYPE<T> >::iterator p;
779 int pInt = 0;
780 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
781 // for (auto& p : _particles) {
782 if (p->getActive()) {
783 frc = p->getForce();
784 vel = p->getVel();
785 p->resetForce();
786 for (auto f : _forces) {
787 f->applyForce(p, pInt, *this);
788 }
789 for (int i = 0; i < 3; i++) {
790 vel[i] += .5 * (p->getForce()[i] + frc[i]) / p->getMass() * dT;
791 }
792 p->setVel(vel);
793 }
794 }
795 }
796
797 /*
798 template<typename T, template<typename U> class PARTICLETYPE>
799 void ParticleSystem3D<T, PARTICLETYPE>::rungeKutta4(T dT)
800 {
801 rungeKutta4_1(dT);
802 rungeKutta4_2(dT);
803 rungeKutta4_3(dT);
804 rungeKutta4_4(dT);
805 }
806
807
808 template<typename T, template<typename U> class PARTICLETYPE>
809 void ParticleSystem3D<T, PARTICLETYPE>::rungeKutta4_1(T dT)
810 {
811 typename std::deque<PARTICLETYPE<T> >::iterator p;
812 int pInt = 0;
813 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
814 if (p->getActive()) {
815 p->resetForce();
816 for (auto f : _forces) {
817 f->applyForce(p, pInt, *this);
818 }
819
820 std::vector<T> k1 = p->getForce();
821 p->setForce(k1, 3);
822 std::vector<T> storeVel = p->getVel();
823 p->setVel(storeVel, 3);
824 std::vector<T> storePos = p->getPos();
825 p->setVel(storePos, 2);
826
827 std::vector<T> vel(3, T()), pos(3, T());
828 for (int i = 0; i < 3; i++) {
829 vel[i] = p->getVel(3)[i] + dT / 2. / p->getMass() * k1[i];
830 pos[i] = p->getVel(2)[i] + dT / 2. * vel[i];
831 }
832 p->setVel(vel);
833 p->setPos(pos);
834 }
835 }
836 }
837
838
839 template<typename T, template<typename U> class PARTICLETYPE>
840 void ParticleSystem3D<T, PARTICLETYPE>::rungeKutta4_2(T dT)
841 {
842 typename std::deque<PARTICLETYPE<T> >::iterator p;
843 int pInt = 0;
844 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
845 if (p->getActive()) {
846 p->resetForce();
847 for (auto f : _forces) {
848 f->applyForce(p, pInt, *this);
849 }
850 std::vector<T> k2 = p->getForce();
851 p->setForce(k2, 2);
852
853 std::vector<T> vel(3, T()), pos(3, T());
854 for (int i = 0; i < 3; i++) {
855 vel[i] = p->getVel(3)[i] + dT / 2. / p->getMass() * k2[i];
856 pos[i] = p->getVel(2)[i] + dT / 2. * vel[i];
857 }
858 p->setVel(vel);
859 p->setPos(pos);
860 }
861 }
862 }
863
864 template<typename T, template<typename U> class PARTICLETYPE>
865 void ParticleSystem3D<T, PARTICLETYPE>::rungeKutta4_3(T dT)
866 {
867 typename std::deque<PARTICLETYPE<T> >::iterator p;
868 int pInt = 0;
869 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
870 if (p->getActive()) {
871 p->resetForce();
872 for (auto f : _forces) {
873 f->applyForce(p, pInt, *this);
874 }
875 std::vector<T> k3 = p->getForce();
876 p->setForce(k3, 1);
877 std::vector<T> vel(3, T()), pos(3, T());
878 for (int i = 0; i < 3; i++) {
879 vel[i] = p->getVel(3)[i] + dT / p->getMass() * k3[i];
880 pos[i] = p->getVel(2)[i] + dT * vel[i];
881 }
882 p->setVel(vel);
883 p->setPos(pos);
884 }
885 }
886 }
887
888 template<typename T, template<typename U> class PARTICLETYPE>
889 void ParticleSystem3D<T, PARTICLETYPE>::rungeKutta4_4(T dT)
890 {
891 typename std::deque<PARTICLETYPE<T> >::iterator p;
892 int pInt = 0;
893 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
894 if (p->getActive()) {
895 p->resetForce();
896 for (auto f : _forces) {
897 f->applyForce(p, pInt, *this);
898 }
899 std::vector<T> k4 = p->getForce();
900
901 std::vector<T> vel(3, T()), pos(3, T());
902 for (int i = 0; i < 3; i++) {
903 vel[i] = p->getVel(3)[i]
904 + dT / 6. / p->getMass()
905 * (p->getForce(3)[i] + 2 * p->getForce(2)[i]
906 + 2 * p->getForce(1)[i] + p->getForce()[i]);
907 pos[i] = p->getVel(2)[i] + dT * vel[i];
908 }
909 p->setVel(vel);
910 p->setPos(pos);
911 }
912 }
913 }
914 */
915
916 template<typename T, template<typename U> class PARTICLETYPE>
917 std::deque<PARTICLETYPE<T>*> ParticleSystem3D<T, PARTICLETYPE>::
918 getParticlesPointer()
919 {
920 std::deque<PARTICLETYPE<T>*> pointerParticles;
921 for (int i = 0; i < (int) _particles.size(); i++) {
922 pointerParticles.push_back(&_particles[i]);
923 }
924 return pointerParticles;
925 }
926
927 template<typename T, template<typename U> class PARTICLETYPE>
928 std::list<std::shared_ptr<Force3D<T, PARTICLETYPE> > >
933
934 template<typename T, template<typename U> class PARTICLETYPE>
935 std::deque<PARTICLETYPE<T>*> ParticleSystem3D<T, PARTICLETYPE>::
936 getAllParticlesPointer()
937 {
938 std::deque<PARTICLETYPE<T>*> pointerParticles;
939 int i = 0;
940 for (; i < (int) _particles.size(); i++) {
941 pointerParticles.push_back(&_particles[i]);
942 }
943 for (; i < (int) (_particles.size() + _shadowParticles.size()); i++) {
944 pointerParticles.push_back(&_shadowParticles[i - _particles.size()]);
945 }
946 return pointerParticles;
947 }
948
949 template<typename T, template<typename U> class PARTICLETYPE>
950 std::deque<PARTICLETYPE<T>*> ParticleSystem3D<T, PARTICLETYPE>::
951 getShadowParticlesPointer()
952 {
953 std::deque<PARTICLETYPE<T>*> pointerParticles;
954 for (int i = 0; i < (int) (_shadowParticles.size()); i++) {
955 pointerParticles.push_back(&_shadowParticles[i]);
956 }
957 return pointerParticles;
958 }
959
960 template<typename T, template<typename U> class PARTICLETYPE>
961 void ParticleSystem3D<T, PARTICLETYPE>::saveToFile(std::string fullName)
962 {
963 std::ofstream fout(fullName.c_str(), std::ios::app);
964 if (!fout) {
965 clout << "Error: could not open " << fullName << std::endl;
966 }
967 std::vector<T> serPar(PARTICLETYPE<T>::serialPartSize, 0);
968 for (auto& p : _particles) {
969 p.serialize(&serPar[0]);
970 typename std::vector<T>::iterator it = serPar.begin();
971 for (; it != serPar.end(); ++it) {
972 fout << *it << " ";
973 }
974 fout << std::endl;
975 //fout << p.getPos()[0] << " " << p.getPos()[1] << " " << p.getPos()[2] << " " << p.getVel()[0] << " " << p.getVel()[1] << " "<< p.getVel()[2] << std::endl;
976 }
977 fout.close();
978 }
979
980 //template<typename T, template<typename U> class PARTICLETYPE>
981 //template<template<typename V> class DESCRIPTOR>
982 //void ParticleSystem3D<T, PARTICLETYPE>::particleOnFluid(
983 // BlockLattice<T, DESCRIPTOR>& bLattice, Cuboid3D<T>& cuboid,
984 // int overlap, T eps, BlockGeometry<T,3>& bGeometry)
985 //{
986 // T rad = 0;
987 // std::vector<T> vel(3, T());
988 // T minT[3] = {0}, maxT[3] = {0}, physR[3] = {0};
989 // int min[3] = {0}, max[3] = {0};
990 // //cout << "OVERLAP: " << overlap << std::endl;
991 // for (int i = 0; i < sizeInclShadow(); ++i) {
992 // rad = operator[](i).getRad();
993 // minT[0] = operator[](i).getPos()[0] - rad * eps;
994 // minT[1] = operator[](i).getPos()[1] - rad * eps;
995 // minT[2] = operator[](i).getPos()[2] - rad * eps;
996 // maxT[0] = operator[](i).getPos()[0] + rad * eps;
997 // maxT[1] = operator[](i).getPos()[1] + rad * eps;
998 // maxT[2] = operator[](i).getPos()[2] + rad * eps;
999 // cuboid.getLatticeR(min, minT);
1000 // cuboid.getLatticeR(max, maxT);
1001 // max[0]++;
1002 // max[1]++;
1003 // max[2]++;
1004 // T porosity = 0;
1005 // T dist = 0;
1006 // for (int iX = min[0]; iX < max[0]; ++iX) {
1007 // for (int iY = min[1]; iY < max[1]; ++iY) {
1008 // for (int iZ = min[2]; iZ < max[2]; ++iZ) {
1009 // cuboid.getPhysR(physR, {iX, iY, iZ});
1010 // dist = util::pow(physR[0] - operator[](i).getPos()[0], 2)
1011 // + util::pow(physR[1] - operator[](i).getPos()[1], 2)
1012 // + util::pow(physR[2] - operator[](i).getPos()[2], 2);
1013 // if (dist < rad * rad * eps * eps) {
1014 // vel = operator[](i).getVel();
1015 // vel[0] = _converter.latticeVelocity(vel[0]);
1016 // vel[1] = _converter.latticeVelocity(vel[1]);
1017 // vel[2] = _converter.latticeVelocity(vel[2]);
1018 // if (dist < rad * rad) {
1019 // porosity = 0;
1020 // bLattice.get(iX, iY, iZ).defineField<descriptors::POROSITY>(
1021 // &porosity);
1022 // bLattice.get(iX, iY, iZ).defineField<descriptors::LOCAL_DRAG>(
1023 // &vel[0]);
1024 // } else {
1025 // T d = util::sqrt(dist) - rad;
1026 // porosity = 1.
1027 // - util::pow(util::cos(M_PI * d / (2. * (eps - 1.) * rad)), 2);
1028 // bLattice.get(iX, iY, iZ).defineField<descriptors::POROSITY>(
1029 // &porosity);
1030 // bLattice.get(iX, iY, iZ).defineField<descriptors::LOCAL_DRAG>(
1031 // &vel[0]);
1032 // }
1033 // }
1034 // }
1035 // }
1036 // }
1037 // }
1038 //}
1039 //
1040 //template<typename T, template<typename U> class PARTICLETYPE>
1041 //template<template<typename V> class DESCRIPTOR>
1042 //void ParticleSystem3D<T, PARTICLETYPE>::resetFluid(
1043 // BlockLattice<T, DESCRIPTOR>& bLattice, Cuboid3D<T>& cuboid,
1044 // int overlap)
1045 //{
1046 // T rad = 0, vel[3] = {0};
1047 // T minT[3] = {0}, maxT[3] = {0}, physR[3] = {0};
1048 // int min[3] = {0}, max[3] = {0};
1049 // for (int i = 0; i < sizeInclShadow(); ++i) {
1050 // rad = operator[](i).getRad();
1051 // minT[0] = operator[](i).getPos()[0] - rad * 1.2;
1052 // minT[1] = operator[](i).getPos()[1] - rad * 1.2;
1053 // minT[2] = operator[](i).getPos()[2] - rad * 1.2;
1054 // maxT[0] = operator[](i).getPos()[0] + rad * 1.2;
1055 // maxT[1] = operator[](i).getPos()[1] + rad * 1.2;
1056 // maxT[2] = operator[](i).getPos()[2] + rad * 1.2;
1057 // cuboid.getLatticeR(min, minT);
1058 // cuboid.getLatticeR(max, maxT);
1059 // max[0]++;
1060 // max[1]++;
1061 // max[2]++;
1062 // T porosity = 1;
1063 // T dist = 0;
1064 // for (int iX = min[0]; iX < max[0]; ++iX) {
1065 // for (int iY = min[1]; iY < max[1]; ++iY) {
1066 // for (int iZ = min[2]; iZ < max[2]; ++iZ) {
1067 // bLattice.get(iX, iY, iZ).defineField<descriptors::POROSITY>(
1068 // &porosity);
1069 // bLattice.get(iX, iY, iZ).defineField<descriptors::LOCAL_DRAG>(
1070 // vel);
1071 // // }
1072 // }
1073 // }
1074 // }
1075 // }
1076 //}
1077
1078
1079 template<>
1080 void ParticleSystem3D<double, MagneticParticle3D>::resetMag()
1081 {
1082 typename std::deque<MagneticParticle3D<double> >::iterator p;
1083 int pInt = 0;
1084 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1085
1086 if (p->getActive()) {
1087
1088 p->resetForce();
1089 p->resetTorque();
1090 }
1091 }
1092 }
1093
1094 // multiple collision models
1095 template<>
1096 void ParticleSystem3D<double, MagneticParticle3D>::resetMag(std::set<int> sActivityOfParticle)
1097 {
1098 typename std::deque<MagneticParticle3D<double> >::iterator p;
1099 int pInt = 0;
1100 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1101
1102 if (p->getSActivity() == 3) {
1103 continue;
1104 }
1105
1106 if (p->getActive()) {
1107 bool b = false;
1108 for (auto sA : sActivityOfParticle) {
1109 if (p->getSActivity() == sA) {
1110 b = true;
1111 break;
1112 }
1113 }
1114 if (b == false) {
1115 continue;
1116 }
1117 p->resetForce();
1118 p->resetTorque();
1119 }
1120 }
1121 }
1122
1123
1124 template<>
1125 void ParticleSystem3D<double, MagneticParticle3D>::computeForce()
1126 {
1127 typename std::deque<MagneticParticle3D<double> >::iterator p;
1128 int pInt = 0;
1129 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1130 if (p->getActive()) {
1131
1132 for (auto f : _forces) {
1133 f->applyForce(p, pInt, *this);
1134 }
1135 }
1136 }
1137 }
1138
1139 // multiple collision models
1140 template<>
1141 void ParticleSystem3D<double, MagneticParticle3D>::computeForce(std::set<int> sActivityOfParticle)
1142 {
1143 typename std::deque<MagneticParticle3D<double> >::iterator p;
1144 int pInt = 0;
1145 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1146
1147 if (p->getActive()) {
1148
1149 if (p->getSActivity() == 3) {
1150 continue;
1151 }
1152
1153 bool b = false;
1154 for (auto sA : sActivityOfParticle) {
1155 if (p->getSActivity() == sA) {
1156 b = true;
1157 break;
1158 }
1159 }
1160 if (b == false) {
1161 continue;
1162 }
1163
1164 for (auto f : _forces) {
1165 // f->applyForce(p, p->getID(), *this);
1166 f->applyForce(p, pInt, *this);
1167 }
1168 }
1169 }
1170 }
1171
1172 /* Original MagDM damping intern
1173 template<>
1174 void ParticleSystem3D<double, MagneticParticle3D>::integrateTorqueMag(double dT)
1175 {
1176 for (auto& p : _particles) {
1177 Vector<double, 3> deltaAngle;
1178 double angle;
1179 double epsilon = std::numeric_limits<double>::epsilon();
1180 double damping = util::pow((1. - p.getADamping()), dT);
1181 for (int i = 0; i < 3; i++) {
1182 p.getAVel()[i] += (5. * (p.getTorque()[i]) * dT) / (2. * p.getMass() * util::pow(p.getRad(), 2));
1183 p.getAVel()[i] *= damping;
1184 deltaAngle[i] = p.getAVel()[i] * dT;
1185
1186 angle = norm(deltaAngle);
1187 if (angle > epsilon) {
1188 std::vector<double> null(3, double());
1189
1190 RotationRoundAxis3D<double, double> rotRAxis(null, util::fromVector3(deltaAngle), angle);
1191 double input[3] = {p.getMoment()[0], p.getMoment()[1], p.getMoment()[2]};
1192 Vector<double, 3> in(input);
1193 double output[3] = {double(), double(), double()};
1194 rotRAxis(output, input);
1195 Vector<double, 3> out(output);
1196 // renormalize output
1197 if (out.norm() > epsilon) {
1198 out = (1. / out.norm()) * out;
1199 }
1200
1201 p.getMoment()[0] = out[0];
1202 p.getMoment()[1] = out[1];
1203 p.getMoment()[2] = out[2];
1204 }
1205 }
1206 }
1207 }
1208 */
1209
1210 template<>
1211 void ParticleSystem3D<double, MagneticParticle3D>::integrateTorqueMag(double dT)
1212 {
1213 for (auto& p : _particles) {
1214
1215 Vector<double, 3> deltaAngle;
1216 double angle;
1217 double epsilon = std::numeric_limits<double>::epsilon();
1218 for (int i = 0; i < 3; i++) {
1219 p.getAVel()[i] += (5. * (p.getTorque()[i]) * dT) / (2. * p.getMass() * util::pow(p.getRad(), 2));
1220 deltaAngle[i] = p.getAVel()[i] * dT;
1221 angle = norm(deltaAngle);
1222
1223 if (angle > epsilon) {
1224 std::vector<double> null(3, double());
1225 RotationRoundAxis3D<double, double> rotRAxis(null, util::fromVector3(deltaAngle), angle);
1226 double input[3] = {p.getMoment()[0], p.getMoment()[1], p.getMoment()[2]};
1227 Vector<double, 3> in(input);
1228 double output[3] = {double(), double(), double()};
1229 rotRAxis(output, input);
1230 Vector<double, 3> out(output);
1231 // renormalize output
1232 if (norm(out) > epsilon) {
1233 out = normalize(out);
1234 }
1235
1236 p.getMoment()[0] = out[0];
1237 p.getMoment()[1] = out[1];
1238 p.getMoment()[2] = out[2];
1239 }
1240 }
1241 }
1242 }
1243
1244 // multiple collision models
1245 template<>
1246 void ParticleSystem3D<double, MagneticParticle3D>::integrateTorqueMag(double dT, std::set<int> sActivityOfParticle)
1247 {
1248
1249 for (auto& p : _particles) {
1250
1251 if (p.getSActivity() == 3) {
1252 continue;
1253 }
1254
1255 bool b = false;
1256 for (auto sA : sActivityOfParticle) {
1257 if (p.getSActivity() == sA) {
1258 b = true;
1259 break;
1260 }
1261 }
1262 if (b == false) {
1263 continue;
1264 }
1265
1266 Vector<double, 3> deltaAngle;
1267 double angle;
1268 double epsilon = std::numeric_limits<double>::epsilon();
1269 for (int i = 0; i < 3; i++) {
1270 p.getAVel()[i] += (5. * (p.getTorque()[i]) * dT) / (2. * p.getMass() * util::pow(p.getRad(), 2));
1271 deltaAngle[i] = p.getAVel()[i] * dT;
1272 angle = norm(deltaAngle);
1273
1274 if (angle > epsilon) {
1275 std::vector<double> null(3, double());
1276 RotationRoundAxis3D<double, double> rotRAxis(null, util::fromVector3(deltaAngle), angle);
1277 double input[3] = {p.getMoment()[0], p.getMoment()[1], p.getMoment()[2]};
1278 Vector<double, 3> in(input);
1279 double output[3] = {double(), double(), double()};
1280 rotRAxis(output, input);
1281 Vector<double, 3> out(output);
1282 // renormalize output
1283 if (norm(out) > epsilon) {
1284 out = normalize(out);
1285 }
1286
1287 p.getMoment()[0] = out[0];
1288 p.getMoment()[1] = out[1];
1289 p.getMoment()[2] = out[2];
1290 }
1291 }
1292 }
1293 }
1294
1295 template<typename T, template<typename U> class PARTICLETYPE>
1296 void ParticleSystem3D<T, PARTICLETYPE>::setStoredValues()
1297 {
1298
1299 typename std::deque<PARTICLETYPE<T> >::iterator p;
1300 int pInt = 0;
1301 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1302
1303 p->setStoredPos(p->getPos());
1304 // p->setStoredVel(p->getVel());
1305 // p->setStoreForce(p->getForce());
1306 }
1307 }
1308
1309 template<>
1310 void ParticleSystem3D<double, MagneticParticle3D>::setOverlapZero()
1311 {
1312
1313 typename std::deque<MagneticParticle3D<double> >::iterator p;
1314 int pInt = 0;
1315 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1316
1317 std::vector<std::pair<size_t, double>> ret_matches;
1318 // kind of contactDetection has to be chosen in application
1319 getContactDetection()->getMatches(pInt, ret_matches);
1320
1321 MagneticParticle3D<double>* p2 = NULL;
1322
1323 // iterator walks through number of neighbored particles = ret_matches
1324 for (const auto& it : ret_matches) {
1325
1326 if (!util::nearZero(it.second)) {
1327 p2 = &(_particles.at(it.first)); //p2 = &pSys[it.first];
1328
1329 if ((p2->getRad() + p->getRad()) > (util::sqrt(it.second))) {
1330
1331 // overlap
1332 double overlap = (p2->getRad() + p->getRad()) - util::sqrt(it.second);
1333
1334 //conVec: vector from particle1 to particle2
1335 Vector<double, 3> conVec(0., 0., 0.) ;
1336 for (int i = 0; i <= 2; i++) {
1337
1338 conVec[i] = p2->getPos()[i] - p->getPos()[i];
1339 }
1340 Vector<double, 3> conVecNormalized(conVec) ;
1341 normalize(conVecNormalized) ;
1342
1343 double dpos[3] = {double(0), double(0), double(0) } ;
1344
1345 // Both particles are not deposited (sActivity = 3)
1346 if ((p->getSActivity() != 3) && (p2->getSActivity() != 3)) {
1347
1348 for (int i = 0; i <= 2; i++) {
1349
1350 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1351 p->getPos()[i] -= 1.* dpos[i];
1352 p2->getPos()[i] += 1.* dpos[i];
1353 }
1354 if ((p->getSActivity() == 2) || (p->getSActivity() == 2)) {
1355 p->setSActivity(2);
1356 p2->setSActivity(2);
1357 }
1358 continue;
1359 }
1360
1361 // Particle 1 is deposited (sActivity = 3) and Particle 2 is not
1362 if ((p->getSActivity() != 3) && (p2->getSActivity() == 3)) {
1363
1364 for (int i = 0; i <= 2; i++) {
1365
1366 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1367 p->getPos()[i] -= 2. * dpos[i];
1368 }
1369 p->setSActivity(2) ;
1370 }
1371
1372 // Particle 2 is deposited (sActivity = 3) and Particle 1 is not
1373 if ((p->getSActivity() == 3) && (p2->getSActivity() != 3)) {
1374
1375 for (int i = 0; i <= 2; i++) {
1376
1377 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1378 p2->getPos()[i] += 2. * dpos[i];
1379 }
1380 p2->setSActivity(2) ;
1381 }
1382
1383 // Both particles are not deposited (sActivity = 3)
1384 if ((p->getSActivity() == 3) && (p2->getSActivity() == 3)) {
1385
1386 for (int i = 0; i <= 2; i++) {
1387
1388 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1389 p->getPos()[i] -= dpos[i];
1390 p2->getPos()[i] += dpos[i];
1391 }
1392 }
1393 }
1394 }
1395 }
1396 }
1397 }
1398
1399 template<>
1400 void ParticleSystem3D<double, MagneticParticle3D>::setOverlapZeroForCombinationWithMechContactForce()
1401 {
1402
1403 typename std::deque<MagneticParticle3D<double> >::iterator p;
1404 int pInt = 0;
1405
1406 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1407
1408 if (p->getSActivity() > 1) {
1409 continue;
1410 }
1411
1412 std::vector<std::pair<size_t, double>> ret_matches;
1413 // kind of contactDetection has to be chosen in application
1414 getContactDetection()->getMatches(pInt, ret_matches);
1415
1416 MagneticParticle3D<double>* p2 = NULL;
1417 // iterator walks through number of neighbored particles = ret_matches
1418 for (const auto& it : ret_matches) {
1419 if (!util::nearZero(it.second)) {
1420 p2 = &(_particles.at(it.first)); //p2 = &pSys[it.first];
1421
1422 if ((p2->getRad() + p->getRad()) > (util::sqrt(it.second))) {
1423 // overlap
1424 double overlap = (p2->getRad() + p->getRad()) - util::sqrt(it.second);
1425
1426 //conVec: vector from particle1 to particle2
1427 Vector<double, 3> conVec(0., 0., 0.) ;
1428 for (int i = 0; i <= 2; i++) {
1429
1430 conVec[i] = p2->getPos()[i] - p->getPos()[i];
1431 }
1432 Vector<double, 3> conVecNormalized(conVec) ;
1433 normalize(conVecNormalized) ;
1434
1435 double dpos[3] = {double(0), double(0), double(0) } ;
1436
1437 // Both particles are not deposited
1438 if (overlap > p2->getRad() + p->getRad()) {
1439
1440 if (p2->getSActivity() == 1) {
1441
1442 for (int i = 0; i <= 2; i++) {
1443
1444 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1445 p->getPos()[i] -= dpos[i] * 1. ;
1446 p2->getPos()[i] += dpos[i] * 1. ;
1447 }
1448 continue;
1449 }
1450 }
1451
1452 // Particle 2 is out of the sphere of influence of setOverlapZeroForCombinationWithMechContactForce()
1453 // Particle 1 is transferred to the influence of the mechanic contact force
1454 else {
1455 for (int i = 0; i <= 2; i++) {
1456
1457 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1458 p->getPos()[i] -= 2 * dpos[i];
1459 }
1460 p->setSActivity(2);
1461 continue;
1462 }
1463 }
1464 }
1465 }
1466 }
1467 }
1468
1469 template<typename T, template<typename U> class PARTICLETYPE>
1470 void ParticleSystem3D<T, PARTICLETYPE>::getMinDistParticle (std::vector<std::pair<size_t, T>> ret_matches)
1471 {
1472 std::sort(ret_matches.begin(), ret_matches.end(), getMinDistPartObj);
1473 }
1474
1475 template<>
1476 void ParticleSystem3D<double, MagneticParticle3D>::partialElasticImpact(double restitutionCoeff)
1477 {
1478 typename std::deque<MagneticParticle3D<double> >::iterator p;
1479 int pInt = 0;
1480
1481 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1482
1483 std::vector<std::pair<size_t, double>> ret_matches;
1484 // kind of contactDetection has to be chosen in application
1485 getContactDetection()->getMatches(pInt, ret_matches);
1486
1487 MagneticParticle3D<double>* p2 = NULL;
1488
1489 // iterator walks through number of neighbored particles = ret_matches
1490 for (const auto& it : ret_matches) {
1491
1492 if (!util::nearZero(it.second)) {
1493 p2 = &(_particles.at(it.first));
1494
1495 if ((p2->getRad() + p->getRad()) > (util::sqrt(it.second))) {
1496
1497 // overlap
1498 double overlap = (p2->getRad() + p->getRad()) - util::sqrt(it.second);
1499
1500 //conVec: vector from particle1 to particle2
1501 Vector<double, 3> conVec(0., 0., 0.) ;
1502 for (int i = 0; i <= 2; i++) {
1503
1504 conVec[i] = p2->getPos()[i] - p->getPos()[i];
1505 }
1506 Vector<double, 3> conVecNormalized(conVec) ;
1507 normalize(conVecNormalized) ;
1508
1509 // Particle velocities before collision
1510 Vector<double, 3> vel1bc = {p->getVel()} ;
1511 Vector<double, 3> vel2bc = {p2->getVel()} ;
1512
1513 // Normal and tangential particle velocities before collision
1514 Vector<double, 3> velN1(0., 0., 0.) ;
1515 Vector<double, 3> velT1(0., 0., 0.) ;
1516 Vector<double, 3> velN2(0., 0., 0.) ;
1517 Vector<double, 3> velT2(0., 0., 0.) ;
1518
1519 // Particle velocities after collision
1520 Vector<double, 3> vel1ac(0., 0., 0.) ;
1521 Vector<double, 3> vel2ac(0., 0., 0.) ;
1522
1523 // Restitution coefficient
1524 double Cr = restitutionCoeff;
1525
1526 velN1 = (conVecNormalized * vel1bc) * conVecNormalized;
1527 velN2 = (conVecNormalized * vel2bc) * conVecNormalized;
1528 velT1 = vel1bc - velN1;
1529 velT2 = vel2bc - velN2;
1530 vel1ac = (Cr * p2->getMass() * ( velN2 - velN1) + (p2->getMass() * velN2) + (p->getMass() * velN1)) * (1. / (p->getMass() + p2->getMass())) ;
1531 vel2ac = (Cr * p->getMass() * ( velN1 - velN2) + (p2->getMass() * velN2) + (p->getMass() * velN1)) * (1. / (p->getMass() + p2->getMass())) ;
1532
1533 double dpos[3] = {double(0), double(0), double(0) } ;
1534
1535 // Both particles are not deposited (sActivity = 3)
1536 if ((p->getSActivity() != 3) && (p2->getSActivity() != 3)) {
1537
1538 for (int i = 0; i <= 2; i++) {
1539
1540 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1541 p->getPos()[i] -= dpos[i] * 1. ;
1542 p->getVel()[i] = vel1ac[i] + velT1[i] ;
1543 p2->getPos()[i] += dpos[i] * 1. ;
1544 p2->getVel()[i] = vel2ac[i] + velT2[i] ;
1545 }
1546 if ((p->getSActivity() == 2) || (p->getSActivity() == 2)) {
1547 p->setSActivity(2);
1548 p2->setSActivity(2);
1549 }
1550 continue;
1551 }
1552
1553 // Particle 1 is deposited (sActivity = 3) and Particle 2 is not
1554 if ((p->getSActivity() != 3) && (p2->getSActivity() == 3)) {
1555
1556 for (int i = 0; i <= 2; i++) {
1557
1558 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1559 p->getPos()[i] -= 2. * dpos[i];
1560 p->getVel()[i] = -1. * p->getVel()[i];
1561 }
1562 p->setSActivity(2) ;
1563 continue;
1564 }
1565
1566 // Particle 2 is deposited (sActivity = 3) and Particle 1 is not
1567 if ((p->getSActivity() == 3) && (p2->getSActivity() != 3)) {
1568
1569 for (int i = 0; i <= 2; i++) {
1570
1571 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1572 p2->getPos()[i] += 2. * dpos[i];
1573 p2->getVel()[i] = -1. * p2->getVel()[i];
1574 }
1575 p2->setSActivity(2) ;
1576 continue;
1577 }
1578
1579 // Both particles are deposited (sActivity = 3)
1580 if ((p->getSActivity() == 3) && (p2->getSActivity() == 3)) {
1581
1582 for (int i = 0; i <= 2; i++) {
1583
1584 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1585 p->getPos()[i] -= dpos[i];
1586 p2->getPos()[i] += dpos[i];
1587 }
1588 continue;
1589 }
1590 }
1591 }
1592 }
1593 }
1594 }
1595
1596 template<>
1597 void ParticleSystem3D<double, MagneticParticle3D>::partialElasticImpactV2(double restitutionCoeff)
1598 {
1599 typename std::deque<MagneticParticle3D<double> >::iterator p;
1600 int pInt = 0;
1601
1602 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1603
1604 std::vector<std::pair<size_t, double>> ret_matches;
1605 // kind of contactDetection has to be chosen in application
1606 getContactDetection()->getMatches(pInt, ret_matches);
1607
1608 MagneticParticle3D<double>* p2 = NULL;
1609
1610 // iterator walks through number of neighbored particles = ret_matches
1611 if (ret_matches.size() == 1) {
1612 continue;
1613 }
1614 double minDist = 0. ;
1615 int xPInt = 0 ;
1616
1617 for (auto a : ret_matches) {
1618 if (util::nearZero(a.second)) {
1619 continue;
1620 }
1621 if (minDist == 0) {
1622 minDist = a.second;
1623 xPInt = a.first;
1624 }
1625 else {
1626 if (a.second < minDist) {
1627 minDist = a.second;
1628 xPInt = a.first;
1629 }
1630 continue;
1631 }
1632 }
1633
1634 p2 = &(_particles.at(xPInt));
1635
1636 if ((p2->getRad() + p->getRad()) > (util::sqrt(minDist))) {
1637
1638 // overlap
1639 double overlap = (p2->getRad() + p->getRad()) - util::sqrt(minDist);
1640
1641 //conVec: vector from particle 1 to particle 2
1642 Vector<double, 3> conVec(0., 0., 0.) ;
1643 for (int i = 0; i <= 2; i++) {
1644
1645 conVec[i] = p2->getPos()[i] - p->getPos()[i];
1646 }
1647 Vector<double, 3> conVecNormalized(conVec) ;
1648 normalize(conVecNormalized) ;
1649
1650 // Particle velocities before collision
1651 Vector<double, 3> vel1bc = {p->getVel()} ;
1652 Vector<double, 3> vel2bc = {p2->getVel()};
1653
1654 // Normal and tangential particle velocities before collision
1655 Vector<double, 3> velN1(0., 0., 0.) ;
1656 Vector<double, 3> velT1(0., 0., 0.) ;
1657 Vector<double, 3> velN2(0., 0., 0.) ;
1658 Vector<double, 3> velT2(0., 0., 0.) ;
1659
1660 // Particle velocities after collision
1661 Vector<double, 3> vel1ac(0., 0., 0.) ;
1662 Vector<double, 3> vel2ac(0., 0., 0.) ;
1663
1664 // Restitution coeffizient
1665 double Cr = restitutionCoeff;
1666
1667 velN1 = (conVecNormalized * vel1bc) * conVecNormalized;
1668 velN2 = (conVecNormalized * vel2bc) * conVecNormalized;
1669 velT1 = vel1bc - velN1;
1670 velT2 = vel2bc - velN2;
1671 vel1ac = (Cr * p2->getMass() * ( velN2 - velN1) + (p2->getMass() * velN2) + (p->getMass() * velN1)) * (1. / (p->getMass() + p2->getMass())) ;
1672 vel2ac = (Cr * p->getMass() * ( velN1 - velN2) + (p2->getMass() * velN2) + (p->getMass() * velN1)) * (1. / (p->getMass() + p2->getMass())) ;
1673
1674 double dpos[3] = {double(0), double(0), double(0) } ;
1675
1676 // Both particles are not deposited (sActivity = 3)
1677 if ((p->getSActivity() != 3) && (p2->getSActivity() != 3)) {
1678
1679 for (int i = 0; i <= 2; i++) {
1680
1681 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1682 p->getPos()[i] -= dpos[i] * 1. ;
1683 p->getVel()[i] = vel1ac[i] + velT1[i] ;
1684 p2->getPos()[i] += dpos[i] * 1. ;
1685 p2->getVel()[i] = vel2ac[i] + velT2[i] ;
1686 }
1687 if ((p->getSActivity() == 2) || (p->getSActivity() == 2)) {
1688 p->setSActivity(2);
1689 p2->setSActivity(2);
1690 }
1691 continue;
1692 }
1693
1694 // Particle 1 is deposited (sActivity = 3) and Particle 2 is not
1695 if ((p->getSActivity() != 3) && (p2->getSActivity() == 3)) {
1696
1697 for (int i = 0; i <= 2; i++) {
1698
1699 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1700 p->getPos()[i] -= 2. * dpos[i];
1701 p->getVel()[i] = -1. * p->getVel()[i];
1702 }
1703 p->setSActivity(2) ;
1704 continue;
1705 }
1706
1707 // Particle 2 is deposited (sActivity = 3) and Particle 1 is not
1708 if ((p->getSActivity() == 3) && (p2->getSActivity() != 3)) {
1709
1710 for (int i = 0; i <= 2; i++) {
1711
1712 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1713 p2->getPos()[i] += 2. * dpos[i];
1714 p2->getVel()[i] = -1. * p2->getVel()[i];
1715 }
1716 p2->setSActivity(2) ;
1717 continue;
1718 }
1719
1720 // Both particles are not deposited (sActivity = 3)
1721 if ((p->getSActivity() == 3) && (p2->getSActivity() == 3)) {
1722
1723 for (int i = 0; i <= 2; i++) {
1724
1725 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1726 p->getPos()[i] -= dpos[i];
1727 p2->getPos()[i] += dpos[i];
1728 }
1729 continue;
1730 }
1731 }
1732 }
1733 }
1734
1735 template<>
1736 void ParticleSystem3D<double, MagneticParticle3D>::partialElasticImpactForCombinationWithMechContactForce(double restitutionCoeff)
1737 {
1738 typename std::deque<MagneticParticle3D<double> >::iterator p;
1739 int pInt = 0;
1740
1741 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1742
1743 if (p->getSActivity() > 1) {
1744 continue;
1745 }
1746
1747 std::vector<std::pair<size_t, double>> ret_matches;
1748 // kind of contactDetection has to be chosen in application
1749 getContactDetection()->getMatches(pInt, ret_matches);
1750
1751 MagneticParticle3D<double>* p2 = NULL;
1752
1753 // iterator walks through number of neighbored particles = ret_matches
1754 for (const auto& it : ret_matches) {
1755
1756 if (!util::nearZero(it.second)) {
1757 p2 = &(_particles.at(it.first));
1758
1759 if ((p2->getRad() + p->getRad()) > (util::sqrt(it.second))) {
1760
1761 // overlap
1762 double overlap = (p2->getRad() + p->getRad()) - util::sqrt(it.second);
1763
1764 // conVec: vector from particle 1 to particle 2
1765 Vector<double, 3> conVec(0., 0., 0.) ;
1766 for (int i = 0; i <= 2; i++) {
1767
1768 conVec[i] = p2->getPos()[i] - p->getPos()[i];
1769 }
1770
1771 Vector<double, 3> conVecNormalized(conVec) ;
1772 normalize(conVecNormalized) ;
1773
1774 // Particle velocities before collision
1775 Vector<double, 3> vel1bc = {p->getVel()} ;
1776 Vector<double, 3> vel2bc = {p2->getVel()} ;
1777
1778 // Normal and tangential particle velocities before collision
1779 Vector<double, 3> velN1(0., 0., 0.) ;
1780 Vector<double, 3> velT1(0., 0., 0.) ;
1781 Vector<double, 3> velN2(0., 0., 0.) ;
1782 Vector<double, 3> velT2(0., 0., 0.) ;
1783
1784 // Particle velocities after collision
1785 Vector<double, 3> vel1ac(0., 0., 0.) ;
1786 Vector<double, 3> vel2ac(0., 0., 0.) ;
1787
1788 // Restitution coeffizient
1789 double Cr = restitutionCoeff;
1790
1791 velN1 = (conVecNormalized * vel1bc) * conVecNormalized;
1792 velN2 = (conVecNormalized * vel2bc) * conVecNormalized;
1793 velT1 = vel1bc - velN1;
1794 velT2 = vel2bc - velN2;
1795 vel1ac = (Cr * p2->getMass() * ( velN2 - velN1) + (p2->getMass() * velN2) + (p->getMass() * velN1)) * (1. / (p->getMass() + p2->getMass())) ;
1796 vel2ac = (Cr * p->getMass() * ( velN1 - velN2) + (p2->getMass() * velN2) + (p->getMass() * velN1)) * (1. / (p->getMass() + p2->getMass())) ;
1797
1798 double dpos[3] = {double(0), double(0), double(0) } ;
1799
1800 // Both particles are not deposited (sActivity = 3)
1801 if ((p->getSActivity() != 3) && (p2->getSActivity() != 3)) {
1802
1803 for (int i = 0; i <= 2; i++) {
1804
1805 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1806 p->getPos()[i] -= dpos[i] * 1. ;
1807 p->getVel()[i] = vel1ac[i] + velT1[i] ;
1808 p2->getPos()[i] += dpos[i] * 1. ;
1809 p2->getVel()[i] = vel2ac[i] + velT2[i] ;
1810 }
1811 if ((p->getSActivity() == 2) || (p->getSActivity() == 2)) {
1812 p->setSActivity(2);
1813 p2->setSActivity(2);
1814 }
1815 continue;
1816 }
1817
1818 // Particle 1 is deposited (sActivity = 3) and Particle 2 is not
1819 if ((p->getSActivity() != 3) && (p2->getSActivity() == 3)) {
1820
1821 for (int i = 0; i <= 2; i++) {
1822
1823 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1824 p->getPos()[i] -= 2. * dpos[i];
1825 p->getVel()[i] = -1. * p->getVel()[i];
1826 }
1827 p->setSActivity(2) ;
1828 continue;
1829 }
1830
1831 // Particle 2 is deposited (sActivity = 3) and Particle 1 is not
1832 if ((p->getSActivity() == 3) && (p2->getSActivity() != 3)) {
1833
1834 for (int i = 0; i <= 2; i++) {
1835
1836 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1837 p2->getPos()[i] += 2. * dpos[i];
1838 p2->getVel()[i] = -1. * p2->getVel()[i];
1839 }
1840 p2->setSActivity(2) ;
1841 continue;
1842 }
1843
1844 // Both particles are deposited (sActivity = 3)
1845 if ((p->getSActivity() == 3) && (p2->getSActivity() == 3)) {
1846
1847 for (int i = 0; i <= 2; i++) {
1848
1849 dpos[i] = conVecNormalized[i] * 0.5 * overlap ;
1850 p->getPos()[i] -= dpos[i];
1851 p2->getPos()[i] += dpos[i];
1852 }
1853 continue;
1854 }
1855 }
1856 }
1857 }
1858 }
1859 }
1860
1861 template<>
1862 void ParticleSystem3D<double, MagneticParticle3D>::findAgglomerates()
1863 {
1864 typename std::deque<MagneticParticle3D<double> >::iterator p;
1865 int pInt = 0;
1866
1867 for (p = _particles.begin(); p != _particles.end(); ++p, ++pInt) {
1868
1869 auto* p1 = &(*p);
1870 std::vector<std::pair<size_t, double>> ret_matches;
1871
1872 getContactDetection()->getMatches(pInt, ret_matches);
1873
1874 MagneticParticle3D<double>* p2 = NULL;
1875
1876 for (const auto& it : ret_matches) {
1877
1878 if (!util::nearZero(it.second)) {
1879
1880 p2 = &(_particles.at(it.first));
1881
1882 if ((p2->getRad() + p1->getRad()) > (util::sqrt(it.second))) {
1883
1884 // Both particles are non agglomerated
1885 // A new agglomerate is formed
1886 if ((p1->getAggloItr() == _Agglomerates.begin()) && (p2->getAggloItr() == _Agglomerates.begin())) {
1887
1888 std::list<MagneticParticle3D<double>*> aggloList{p1, p2} ;
1889
1890 typename std::list<MagneticParticle3D<double>*>::iterator x1;
1891 typename std::list<MagneticParticle3D<double>*>::iterator x2;
1892
1893 _Agglomerates.push_back(aggloList);
1894 p1->setAggloItr(_Agglomerates.end() - 1);
1895 p2->setAggloItr(_Agglomerates.end() - 1);
1896
1897 for (auto x = _Agglomerates[0].begin(); x != _Agglomerates[0].end(); ++x) {
1898
1899 if (*x == p1) {
1900 x1 = x;
1901 }
1902 if (*x == p2) {
1903 x2 = x;
1904 }
1905
1906 }
1907 _Agglomerates[0].erase(x1);
1908 _Agglomerates[0].erase(x2);
1909
1910 continue;
1911 }
1912
1913 // Particle 2 is already part of an agglomerate and Particle 1 is non agglomerated
1914 // Particle 1 is added to the agglomerate
1915 if ((p1->getAggloItr() == _Agglomerates.begin()) && (p2->getAggloItr() != _Agglomerates.begin())) {
1916
1917 (p2->getAggloItr())->push_back(p1) ;
1918 p1->setAggloItr(p2->getAggloItr()) ;
1919 typename std::list<MagneticParticle3D<double>*>::iterator x1;
1920
1921 for (auto x = _Agglomerates[0].begin(); x != _Agglomerates[0].end(); ++x) {
1922 if (*x == p1) {
1923 x1 = x;
1924 }
1925 }
1926 _Agglomerates[0].erase(x1);
1927
1928 continue;
1929 }
1930
1931 // Particle 1 is already part of an agglomerate and Particle 2 is non agglomerated
1932 // Particle 2 is added to the agglomerate
1933 if ((p1->getAggloItr() != _Agglomerates.begin()) && (p2->getAggloItr() == _Agglomerates.begin())) {
1934
1935 (p1->getAggloItr())->push_back(p2) ;
1936 p2->setAggloItr(p1->getAggloItr()) ;
1937 typename std::list<MagneticParticle3D<double>*>::iterator x2;
1938
1939 for (auto x = _Agglomerates[0].begin(); x != _Agglomerates[0].end(); ++x) {
1940 if (*x == p2) {
1941 x2 = x;
1942 }
1943 }
1944
1945 _Agglomerates[0].erase(x2);
1946
1947 continue;
1948 }
1949
1950 // Both particles are part of diffrent agglomerates
1951 // The two Agglomerates are united
1952 if (((p1->getAggloItr() != _Agglomerates.begin()) && (p2->getAggloItr() != _Agglomerates.begin())) && ( p1->getAggloItr() != p2->getAggloItr())) {
1953
1954 typename std::deque<std::list<MagneticParticle3D<double>*>>::iterator x1;
1955 typename std::deque<std::list<MagneticParticle3D<double>*>>::iterator x2;
1956
1957 if (p1->getAggloItr() <= p2->getAggloItr()) {
1958 x2 = p2->getAggloItr() ;
1959 x1 = p1->getAggloItr() ;
1960
1961 }
1962 else {
1963 x2 = p1->getAggloItr() ;
1964 x1 = p2->getAggloItr() ;
1965 }
1966
1967 x1->splice(x1->end(), *x2) ;
1968
1969 _Agglomerates.erase(x2) ;
1970
1971 for (auto anew = _Agglomerates.begin(); anew != _Agglomerates.end(); ++anew) {
1972
1973 for (auto pnew = anew->begin(); pnew != anew->end(); ++pnew) {
1974
1975 (*pnew)->setAggloItr(anew);
1976 }
1977 }
1978
1979 continue;
1980
1981 }
1982 }
1983 }
1984 }
1985 }
1986 }
1987
1988 template<>
1989 void ParticleSystem3D<double, MagneticParticle3D>::initAggloParticles()
1990 {
1991 typename std::deque<MagneticParticle3D<double> >::iterator p;
1992 static int pInt = 0;
1993
1994 for (p = _particles.begin() + pInt; p != _particles.end(); ++p, ++pInt) {
1995 p->setAggloItr(_Agglomerates.begin());
1996 MagneticParticle3D<double>* pPointer = &(*p);
1997 _Agglomerates.begin()->push_back(pPointer);
1998 }
1999 }
2000
2001
2002} //namespace olb
2003#endif /* PARTICLESYSTEM_3D_HH */
olb::AnalyticalConst
AnalyticalConst: DD -> XD, where XD is defined by value.size()
Definition analyticalF.h:78
olb::BackCouplingModel
Abstact base class for BaseBackCouplingModel.
Definition backCouplingModels.h:40
olb::BlockGeometry
Representation of a block geometry.
Definition blockGeometry.h:56
olb::BlockGeometry::get
std::enable_if_t< sizeof...(L)==D, int > get(L... latticeR) const
Read-only access to a material number.
Definition blockGeometry.h:107
olb::Boundary3D
Prototype for all particle boundaries.
Definition boundary3D.h:43
olb::ContactDetection::generate
virtual ContactDetection< T, PARTICLETYPE > * generate(ParticleSystem3D< T, PARTICLETYPE > &pSys)
Definition contactDetection.h:51
olb::Force3D
Prototype for all particle forces.
Definition force3D.h:43
olb::ForwardCouplingModel
Abstact base class for all the forward-coupling models Its raison d'etre consists of not being temple...
Definition forwardCouplingModels3D.h:41
olb::MagneticParticle3D::setSActivity
void setSActivity(int sActivity)
Definition magneticParticle3D.hh:129
olb::MagneticParticle3D::getAggloItr
std::deque< std::list< MagneticParticle3D< T > * > >::iterator & getAggloItr()
Definition magneticParticle3D.hh:135
olb::MagneticParticle3D::setAggloItr
void setAggloItr(typename std::deque< std::list< MagneticParticle3D< T > * > >::iterator aggloItr)
Definition magneticParticle3D.hh:141
olb::Particle3D::getVel
std::vector< T > & getVel()
Definition particle3D.hh:218
olb::Particle3D::getPos
std::vector< T > & getPos()
Definition particle3D.hh:176
olb::Particle3D::getMass
const T & getMass()
Definition particle3D.hh:423
olb::Particle3D::getRad
const T & getRad()
Definition particle3D.hh:494
olb::ParticleSystem3D::executeBackwardCoupling
bool executeBackwardCoupling(BackCouplingModel< T, PARTICLETYPE > &backwardCoupling, int material, int subSteps=1)
Definition particleSystem3D.hh:409
olb::ParticleSystem3D::partialElasticImpact
void partialElasticImpact(T restitutionCoeff)
Resets existing particle overlaps in the event of a collision and applies the physics of an partial e...
Definition particleSystem3D.h:260
olb::ParticleSystem3D::clearParticles
void clearParticles()
Removes all particles from system.
Definition particleSystem3D.hh:227
olb::ParticleSystem3D::setPosExt
void setPosExt(std::vector< T > physPos, std::vector< T > physExtend)
Set global coordinates and extends of Particlesystem (SI units)
Definition particleSystem3D.hh:114
olb::ParticleSystem3D::setVelToFluidVel
void setVelToFluidVel(SuperLatticeInterpPhysVelocity3D< T, DESCRIPTOR > &)
Set velocity of all particles to fluid velocity.
Definition particleSystem3D.hh:261
olb::ParticleSystem3D::computeForce
void computeForce()
Compute all forces on particles.
Definition particleSystem3D.hh:284
olb::ParticleSystem3D::executeForwardCoupling
bool executeForwardCoupling(ForwardCouplingModel< T, PARTICLETYPE > &forwardCoupling)
Definition particleSystem3D.hh:395
olb::ParticleSystem3D::getForcesPointer
std::list< std::shared_ptr< Force3D< T, PARTICLETYPE > > > getForcesPointer()
returns shared pointer of forces
Definition particleSystem3D.hh:929
olb::ParticleSystem3D::sizeInclShadow
int sizeInclShadow() const
Get number of particles including shadow particles in ParticleSystem.
Definition particleSystem3D.hh:177
olb::ParticleSystem3D::getShadowParticlesPointer
std::deque< PARTICLETYPE< T > * > getShadowParticlesPointer()
returns deque of pointer to all shadow particles contained in a particleSystem3D
Definition particleSystem3D.hh:951
olb::ParticleSystem3D::addParticle
void addParticle(PARTICLETYPE< T > &p)
Add a particle to ParticleSystem.
Definition particleSystem3D.hh:221
olb::ParticleSystem3D::getParticlesPointer
std::deque< PARTICLETYPE< T > * > getParticlesPointer()
returns deque of pointer to particles (not shadow particles) contained in a particleSystem3D
Definition particleSystem3D.hh:918
olb::ParticleSystem3D::setVelToAnalyticalVel
void setVelToAnalyticalVel(AnalyticalConst3D< T, T > &)
Set particle velocity to analytical velocity (e.g. as initial condition.
Definition particleSystem3D.hh:272
olb::ParticleSystem3D::explicitEuler
void explicitEuler(T dT, bool scale=false)
Integration method: explicit Euler if scale = true, velocity is scaled to maximal velocity maximal ve...
Definition particleSystem3D.hh:423
olb::ParticleSystem3D::simulate
virtual void simulate(T deltatime, bool scale=false)
Definition particleSystem3D.hh:338
olb::ParticleSystem3D::setOverlapZero
void setOverlapZero()
Collision models: Todo: enable for parallel mode Resets existing particle overlaps in the event of a ...
Definition particleSystem3D.h:255
olb::ParticleSystem3D::addShadowParticle
void addShadowParticle(PARTICLETYPE< T > &p)
Definition particleSystem3D.hh:233
olb::ParticleSystem3D::getAllParticlesPointer
std::deque< PARTICLETYPE< T > * > getAllParticlesPointer()
returns deque of pointer to all particles (incl.
Definition particleSystem3D.hh:936
olb::ParticleSystem3D::numOfForces
int numOfForces()
Get number of linked forces in ParticleSystem.
Definition particleSystem3D.hh:215
olb::ParticleSystem3D::printDeep
void printDeep(std::string message="")
Definition particleSystem3D.hh:98
olb::ParticleSystem3D::simulateWithTwoWayCoupling_Davide
virtual void simulateWithTwoWayCoupling_Davide(T dT, ForwardCouplingModel< T, PARTICLETYPE > &forwardCoupling, BackCouplingModel< T, PARTICLETYPE > &backCoupling, int material, int subSteps, bool scale)
Definition particleSystem3D.hh:353
olb::ParticleSystem3D::size
int size()
Get number of particles in ParticleSystem.
Definition particleSystem3D.hh:171
olb::ParticleSystem3D::saveToFile
void saveToFile(std::string name)
Save particle positions to file.
Definition particleSystem3D.hh:961
olb::ParticleSystem3D::partialElasticImpactV2
void partialElasticImpactV2(T restitutionCoeff)
Applies the physics of an partial elastic impact while multiple particle overlapping only to the part...
Definition particleSystem3D.h:263
olb::ParticleSystem3D::velocityVerlet1
void velocityVerlet1(T dT)
Integration methods, each need a special template particle.
Definition particleSystem3D.hh:759
olb::ParticleSystem3D::_particles
std::deque< PARTICLETYPE< T > > _particles
Definition particleSystem3D.h:280
olb::ParticleSystem3D::computeBoundary
void computeBoundary()
Compute boundary contact.
Definition particleSystem3D.hh:312
olb::ParticleSystem3D::_forces
std::list< std::shared_ptr< Force3D< T, PARTICLETYPE > > > _forces
Definition particleSystem3D.h:283
olb::ParticleSystem3D::simulateWithTwoWayCoupling_Mathias
virtual void simulateWithTwoWayCoupling_Mathias(T dT, ForwardCouplingModel< T, PARTICLETYPE > &forwardCoupling, BackCouplingModel< T, PARTICLETYPE > &backCoupling, int material, int subSteps, bool scale)
Definition particleSystem3D.hh:344
olb::ParticleSystem3D::addBoundary
void addBoundary(std::shared_ptr< Boundary3D< T, PARTICLETYPE > > pB)
Add a boundary to ParticleSystem.
Definition particleSystem3D.hh:246
olb::ParticleSystem3D::ParticleSystem3D
ParticleSystem3D()=default
Default constructor for ParticleSystem.
olb::ParticleSystem3D::numOfActiveParticles
int numOfActiveParticles()
Get number of active particles in ParticleSystem.
Definition particleSystem3D.hh:183
olb::ParticleSystem3D::setOverlapZeroForCombinationWithMechContactForce
void setOverlapZeroForCombinationWithMechContactForce()
For the combined use of setOverlapZero() and a mechanic contact force.
Definition particleSystem3D.h:257
olb::ParticleSystem3D::computeParticleOperation
void computeParticleOperation()
Compute operations on particles.
Definition particleSystem3D.hh:325
olb::ParticleSystem3D::getPhysPos
const std::vector< T > & getPhysPos()
Get global coordinates and extends of Particlesystem (SI units)
Definition particleSystem3D.hh:122
olb::ParticleSystem3D::addForce
void addForce(std::shared_ptr< Force3D< T, PARTICLETYPE > > pF)
Add a force to ParticleSystem.
Definition particleSystem3D.hh:239
olb::ParticleSystem3D::getMinDistParticle
void getMinDistParticle(std::vector< std::pair< size_t, T > > ret_matches)
Definition particleSystem3D.hh:1470
olb::ParticleSystem3D::countMaterial
int countMaterial(int mat=1)
Get number of particles in vicinity of material number mat.
Definition particleSystem3D.hh:369
olb::ParticleSystem3D::initAggloParticles
void initAggloParticles()
Adds new generated particles to the list of non agglomerated Particles.
Definition particleSystem3D.h:270
olb::ParticleSystem3D::findAgglomerates
void findAgglomerates()
Detects and manages particle agglomerates.
Definition particleSystem3D.h:268
olb::ParticleSystem3D::setContactDetection
void setContactDetection(ContactDetection< T, PARTICLETYPE > &contactDetection)
Set boundary detection algorithm (for future features)
Definition particleSystem3D.hh:107
olb::ParticleSystem3D::operator[]
PARTICLETYPE< T > & operator[](const int i)
Get reference to a particle in the ParticleSystem runs over all particles incl.
Definition particleSystem3D.hh:134
olb::ParticleSystem3D::getContactDetection
ContactDetection< T, PARTICLETYPE > * getContactDetection()
Definition particleSystem3D.hh:92
olb::ParticleSystem3D::setStoredValues
void setStoredValues()
Stores old particle positions - is used in ..._ActExt.
Definition particleSystem3D.hh:1296
olb::ParticleSystem3D::getPhysExtend
const std::vector< T > & getPhysExtend()
Definition particleSystem3D.hh:128
olb::ParticleSystem3D::addParticleOperation
void addParticleOperation(std::shared_ptr< ParticleOperation3D< T, PARTICLETYPE > > pO)
Add an operation to ParticleSystem.
Definition particleSystem3D.hh:253
olb::ParticleSystem3D::partialElasticImpactForCombinationWithMechContactForce
void partialElasticImpactForCombinationWithMechContactForce(T restitutionCoeff)
For the combined use of partialElasticImpact() and a mechanic contact force.
Definition particleSystem3D.h:265
olb::RotationRoundAxis3D
This class saves coordinates of the resulting point after rotation in the output field.
Definition frameChangeF3D.h:356
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
frameChangeF3D.h
This file contains two different classes of functors, in the FIRST part.
olb::util::sqrt
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
olb::util::fromVector3
std::vector< T > fromVector3(const Vector< T, 3 > &vec)
Definition vectorHelpers.h:126
olb::util::max
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:130
olb::util::pow
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112
olb::util::fabs
cpu::simd::Pack< T > fabs(cpu::simd::Pack< T > value)
Definition pack.h:106
olb::util::nearZero
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
olb
Top level namespace for all of OpenLB.
Definition boundaryPostProcessors2D.h:34
olb::norm
constexpr T norm(const ScalarVector< T, D, IMPL > &a)
Euclidean vector norm.
Definition scalarVector.h:66
olb::normalize
constexpr Vector< T, D > normalize(const ScalarVector< T, D, IMPL > &a, T scale=T{1})
Definition vector.h:245
Generated on Tue Mar 5 2024 09:57:52 for OpenLB by doxygen 1.10.0