Initialize at openlb-1-3

1 files changed, 456 insertions, 0 deletions

diff --git a/src/particles/particleSystem3D.h b/src/particles/particleSystem3D.h
new file mode 100755
index 0000000..ad37d86
--- /dev/null
+++ b/src/particles/particleSystem3D.h
@@ -0,0 +1,456 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2015 Thomas Henn
+ *  E-mail contact: info@openlb.net
+ *  The most recent release of OpenLB can be downloaded at
+ *  <http://www.openlb.net/>
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; either version 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public
+ *  License along with this program; if not, write to the Free
+ *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ *  Boston, MA  02110-1301, USA.
+ */
+
+#ifndef PARTICLE_SYSTEM_3D_H
+#define PARTICLE_SYSTEM_3D_H
+
+// Enables particle collision functions
+// CollisionModels: Can be used as an alternative to a mechanic contact force
+// #define CollisionModels
+// CollisionModelsCombindedWithMechContactForce: Can be used in combination with
+// a mechanic contact force
+// #define CollisionModelsCombindedWithMechContactForce
+
+#include <set>
+#include <vector>
+#include <list>
+#include <deque>
+#include "contactDetection/contactDetection.h"
+#include "geometry/superGeometry3D.h"
+#include "core/blockLatticeStructure3D.h"
+#include "forces/force3D.h"
+#include "functors/analytical/analyticalF.h"
+#include "boundaries/boundary3D.h"
+#include "functors/lattice/latticeFrameChangeF3D.h"
+#include "utilities/vectorHelpers.h"
+#include "twoWayCouplings/twoWayCouplings3D.h"
+
+namespace olb {
+
+template<typename T, typename DESCRIPTOR>
+class SuperLatticeInterpPhysVelocity3D;
+
+template<typename T, template<typename U> class PARTICLETYPE>
+class Force3D;
+template<typename T, template<typename U> class PARTICLETYPE>
+class Boundary3D;
+template<typename T, template<typename U> class PARTICLETYPE>
+class ParticleSystem3D;
+template<typename T, template<typename U> class PARTICLETYPE>
+class SuperParticleSystem3D;
+template<typename T, template<typename U> class PARTICLETYPE>
+class SuperParticleSysVtuWriter;
+template<typename T>
+class SuperParticleSysVtuWriterMag;
+
+template<typename T, template<typename U> class PARTICLETYPE>
+class SimulateParticles {
+public:
+  SimulateParticles(ParticleSystem3D<T, PARTICLETYPE>* ps)
+    : _pSys(ps)
+  {
+  }
+  inline void simulate(T dT, bool scale = false)
+  {
+    _pSys->computeForce();
+    _pSys->explicitEuler(dT, scale);
+    //_pSys->rungeKutta4(dT);
+  }
+  inline void simulate(T dT, std::set<int> sActivityOfParticle, bool scale = false) {simulate(dT, scale);}
+  inline void simulateWithTwoWayCoupling_Mathias ( T dT,
+                                    ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling,
+                                    BackCouplingModel<T,PARTICLETYPE>& backCoupling,
+                                    int material, int subSteps = 1, bool scale = false )
+  {
+    for (int iSubStep=1; iSubStep<=subSteps; iSubStep++) {
+      if (! _pSys->executeForwardCoupling(forwardCoupling) ) {
+        std::cout << " on substep " << iSubStep << std::endl;
+        singleton::exit(1);
+      }
+      _pSys->computeForce();
+      _pSys->explicitEuler(dT/(T)(subSteps), scale);
+      //_pSys->rungeKutta4(dT/(T)(subSteps));
+    }
+    _pSys->executeBackwardCoupling(backCoupling, material);
+  }
+  inline void simulateWithTwoWayCoupling_Davide ( T dT,
+                                    ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling,
+                                    BackCouplingModel<T,PARTICLETYPE>& backCoupling,
+                                    int material, int subSteps = 1, bool scale = false )
+  {
+    for (int iSubStep=1; iSubStep<=subSteps; iSubStep++) {
+      if (! _pSys->executeForwardCoupling(forwardCoupling) ) {
+        std::cout << " on substep " << iSubStep << std::endl;
+        singleton::exit(1);
+      }
+      _pSys->executeBackwardCoupling(backCoupling, material, subSteps);
+      _pSys->computeForce();
+      _pSys->explicitEuler(dT/(T)(subSteps), scale);
+      //_pSys->rungeKutta4(dT/(T)(subSteps));
+    }
+  }
+
+private:
+  ParticleSystem3D<T, PARTICLETYPE>* _pSys;
+
+};
+
+
+template<typename T>
+class SimulateParticles<T, RotatingParticle3D> {
+public:
+  SimulateParticles(ParticleSystem3D<T, RotatingParticle3D>* ps)
+    : _pSys(ps)
+  {
+  }
+  inline void simulate(T dT, bool scale = false)
+  {
+    _pSys->computeForce();
+    _pSys->explicitEuler(dT, scale);
+    _pSys->integrateTorque(dT);
+  }
+private:
+  ParticleSystem3D<T, RotatingParticle3D>* _pSys;
+};
+
+template<typename T>
+class SimulateParticles<T, MagneticParticle3D> {
+public:
+  SimulateParticles(ParticleSystem3D<T, MagneticParticle3D>* ps)
+    : _pSys(ps)
+  {
+  }
+
+  inline void simulate(T dT, bool scale = false)
+  {
+    _pSys->resetMag();
+    _pSys->computeForce();
+    _pSys->explicitEuler(dT, scale);
+    _pSys->integrateTorqueMag(dT);
+
+#ifdef CollisionModels
+    _pSys->partialElasticImpact(0.67);
+#endif
+  }
+  inline void simulateWithTwoWayCoupling_Mathias ( T dT,
+                                    ForwardCouplingModel<T,MagneticParticle3D>& forwardCoupling,
+                                    BackCouplingModel<T,MagneticParticle3D>& backCoupling,
+                                    int material, int subSteps = 1, bool scale = false )
+  {
+    for (int iSubStep=1; iSubStep<=subSteps; iSubStep++) {
+      if (! _pSys->executeForwardCoupling(forwardCoupling) ) {
+        std::cout << " on substep " << iSubStep << std::endl;
+        singleton::exit(1);
+      }
+      _pSys->computeForce();
+      _pSys->explicitEuler(dT/(T)(subSteps), scale);
+      //_pSys->rungeKutta4(dT/(T)(subSteps));
+    }
+    _pSys->executeBackwardCoupling(backCoupling, material);
+  }
+  inline void simulateWithTwoWayCoupling_Davide ( T dT,
+                                    ForwardCouplingModel<T,MagneticParticle3D>& forwardCoupling,
+                                    BackCouplingModel<T,MagneticParticle3D>& backCoupling,
+                                    int material, int subSteps = 1, bool scale = false )
+  {
+    for (int iSubStep=1; iSubStep<=subSteps; iSubStep++) {
+      if (! _pSys->executeForwardCoupling(forwardCoupling) ) {
+        std::cout << " on substep " << iSubStep << std::endl;
+        singleton::exit(1);
+      }
+      _pSys->computeForce();
+      _pSys->explicitEuler(dT/(T)(subSteps), scale);
+      //_pSys->rungeKutta4(dT/(T)(subSteps));
+      _pSys->executeBackwardCoupling(backCoupling, material, subSteps);
+    }
+  }
+
+//   multiple collision models
+  inline void simulate(T dT, std::set<int> sActivityOfParticle , bool scale = false)
+  {
+
+    _pSys->resetMag(sActivityOfParticle);
+    _pSys->computeForce(sActivityOfParticle);
+    _pSys->explicitEuler(dT, sActivityOfParticle, scale);
+    _pSys->integrateTorqueMag(dT, sActivityOfParticle);
+
+#ifdef CollisionModels
+    _pSys->partialElasticImpact(0.67);
+#endif
+
+#ifdef CollisionModelsCombindedWithMechContactForce
+    _pSys->partialElasticImpactForCombinationWithMechContactForce(0.67);
+#endif
+  }
+
+private:
+  ParticleSystem3D<T, MagneticParticle3D>* _pSys;
+
+};
+
+
+template<typename T, template<typename U> class PARTICLETYPE>
+class ParticleSystem3D {
+public:
+
+  /// Default constructor for ParticleSystem
+  ParticleSystem3D() = default;
+  /// Constructor for ParticleSystem
+  ParticleSystem3D(int iGeometry, SuperGeometry3D<T>& superGeometry);
+  /// Copy constructor for ParticleSystem
+  ParticleSystem3D(const ParticleSystem3D<T, PARTICLETYPE>& pS);
+  /// Move constructor for ParticleSystem
+  ParticleSystem3D(ParticleSystem3D<T, PARTICLETYPE> && pS);
+  /// Destructor for ParticleSystem
+  virtual ~ParticleSystem3D() {}
+
+  virtual void simulate(T deltatime, bool scale = false);
+  virtual void simulateWithTwoWayCoupling_Mathias ( T dT,
+                                    ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling,
+                                    BackCouplingModel<T,PARTICLETYPE>& backCoupling,
+                                    int material, int subSteps, bool scale );
+  virtual void simulateWithTwoWayCoupling_Davide ( T dT,
+                                    ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling,
+                                    BackCouplingModel<T,PARTICLETYPE>& backCoupling,
+                                    int material, int subSteps, bool scale );
+  // multiple collision models
+  virtual void simulate(T deltatime, std::set<int> sActivityOfParticle, bool scale = false);
+
+  void printDeep(std::string message="");
+
+  /// Get number of particles in ParticleSystem
+  int size();
+  /// Get number of particles including shadow particles in ParticleSystem
+  int sizeInclShadow() const;
+  /// Get number of active particles in ParticleSystem
+  int numOfActiveParticles();
+  /// Get number of linked forces in ParticleSystem
+  int numOfForces();
+  /// Get number of particles in vincinity of material number mat
+  int countMaterial(int mat = 1);
+
+  /// Add a particle to ParticleSystem
+  void addParticle(PARTICLETYPE<T>& p);
+  /// Removes all particles from system
+  void clearParticles();
+  /// Add a force to ParticleSystem
+  void addForce(std::shared_ptr<Force3D<T, PARTICLETYPE> > pF);
+  /// Add a boundary to ParticleSystem
+  void addBoundary(std::shared_ptr<Boundary3D<T, PARTICLETYPE> > pF);
+
+  /// Get reference to a particle in the ParticleSystem
+  /// runs over all particles incl. shadow particles
+  PARTICLETYPE<T>& operator[](const int i);
+  const PARTICLETYPE<T>& operator[](const int i) const;
+
+  /// Set velocity of all particles to fluid velocity
+  template<typename DESCRIPTOR>
+  void setVelToFluidVel(SuperLatticeInterpPhysVelocity3D<T, DESCRIPTOR> &);
+
+  /// Set particle velocity to analytical velocity (e.g. as inital condition
+  void setVelToAnalyticalVel(AnalyticalConst3D<T, T>&);
+
+  /// Set global coordinates and extends of Particlesystem (SI units)
+  void setPosExt(std::vector<T> physPos, std::vector<T> physExtend);
+
+  /// Get global coordinates and extends of Particlesystem (SI units)
+  const std::vector<T>& getPhysPos();
+  const std::vector<T>& getPhysExtend();
+
+  /// Save particle positions to file
+  void saveToFile(std::string name);
+
+  /// Compute all forces on particles
+  void computeForce();
+  // multiple collision models
+  void computeForce(std::set<int> sActivityOfParticle) {computeForce();};
+
+  /// Stores old particle positions - is used in ..._ActExt
+  void setStoredValues();
+
+  /// Sorts the vector of neighbour Particles by increasing distance
+  struct getMinDistPart {
+    bool operator() (std::pair<size_t, T> i, std::pair<size_t, T> j) { return (i.second < j.second);}
+  } getMinDistPartObj;
+  void getMinDistParticle (std::vector<std::pair<size_t, T>> ret_matches);
+
+  /// Compute boundary contact
+  void computeBoundary();
+
+  /// Set boundary detection algorithm (for future features)
+  void setContactDetection(ContactDetection<T, PARTICLETYPE>& contactDetection);
+  ContactDetection<T, PARTICLETYPE>* getContactDetection();
+
+  /// Particle-Fluid interaction for subgrid scale particles
+  //  template<template<typename V> class DESCRIPTOR>
+  //  void particleOnFluid(BlockLatticeStructure3D<T, DESCRIPTOR>& bLattice,
+  //                       Cuboid3D<T>& cuboid, int overlap, T eps,
+  //                       BlockGeometryStructure3D<T>& bGeometry);
+  //  template<template<typename V> class DESCRIPTOR>
+  //  void resetFluid(BlockLatticeStructure3D<T, DESCRIPTOR>& bLattice,
+  //                  Cuboid3D<T>& cuboid, int overlap);
+
+  friend class SuperParticleSystem3D<T, PARTICLETYPE>;
+  friend class SuperParticleSysVtuWriter<T, PARTICLETYPE>;
+  friend class SuperParticleSysVtuWriterMag<T>;
+  friend class SimulateParticles<T, PARTICLETYPE>;
+
+  //std::map<T, int> radiusDistribution();
+
+  /// Integration method: explicit Euler
+  /// if scale = true, velocity is scaled to maximal velocity
+  /// maximal velocity = _superGeometry.getCuboidGeometry().getMaxDeltaR()/dT
+  void explicitEuler(T dT, bool scale = false);
+  // multiple collision models
+  void explicitEuler(T dT, std::set<int> sActivityOfParticle, bool scale = false) {explicitEuler(dT, scale);};
+  bool executeForwardCoupling(ForwardCouplingModel<T,PARTICLETYPE>& forwardCoupling);
+  bool executeBackwardCoupling(BackCouplingModel<T,PARTICLETYPE>& backwardCoupling, int material, int subSteps=1);
+
+  ContactDetection<T, PARTICLETYPE>* getDetection()
+  {
+    return _contactDetection;
+  }
+
+  int getIGeometry()
+  {
+    return _iGeometry;
+  }
+  /// returns deque of pointer to particles (not shadow particles)
+  /// contained in a particleSystem3D
+  std::deque<PARTICLETYPE<T>*> getParticlesPointer();
+  /// returns deque of pointer to all particles (incl. shadow particles)
+  /// contained in a particleSystem3D
+  std::deque<PARTICLETYPE<T>*> getAllParticlesPointer();
+  /// returns deque of pointer to all shadow particles
+  /// contained in a particleSystem3D
+  std::deque<PARTICLETYPE<T>*> getShadowParticlesPointer();
+
+  /// returns deque of particles (no shadow particles)
+  /// contained in a particleSystem3D
+  std::deque<PARTICLETYPE<T>> getParticles()
+  {
+    return _particles;
+  }
+
+  /// returns shared pointer of forces
+  std::list<std::shared_ptr<Force3D<T, PARTICLETYPE> > > getForcesPointer();
+
+  /// Deque of Lists of agglomerated particles
+  std::deque<std::list<PARTICLETYPE<T>*>> _Agglomerates;
+
+protected:
+  void integrateTorque(T dT);
+  void integrateTorqueMag(T dT) {};
+  // multiple collision models
+  void integrateTorqueMag(T dT, std::set<int> sActivityOfParticle) {};
+  void resetMag() {};
+  // multiple collision models
+  void resetMag(std::set<int> sActivityOfParticle) {};
+
+  /// Collision models: Todo: enable for paralle mode
+  /// Resets existing particle overlaps in the event of a collision
+  void setOverlapZero() {};
+  /// For the combined use of setOverlapZero() and a mechanic contact force
+  void setOverlapZeroForCombinationWithMechContactForce() {};
+  /// Resets existing particle overlaps in the event of a collision
+  /// and applies the physics of an partial elastic impact
+  void partialElasticImpact(T restitutionCoeff) {};
+  /// Applies the physics of an partial elastic impact while multiple
+  /// particle overlapping only to the particles with the least separation distance
+  void partialElasticImpactV2(T restitutionCoeff) {};
+  /// For the combined use of partialElasticImpact() and a mechanic contact force
+  void partialElasticImpactForCombinationWithMechContactForce(T restitutionCoeff) {};
+  
+  /// Detects and manages particle agglomerates
+  void findAgglomerates() {};
+  /// Adds new generated particles to the list of non agglomerated Particles
+  void initAggloParticles() {};
+
+  void addShadowParticle(PARTICLETYPE<T>& p);
+
+  mutable OstreamManager clout;
+  int _iGeometry = -1;
+  SuperGeometry3D<T>& _superGeometry;
+  ContactDetection<T, PARTICLETYPE>* _contactDetection;
+  SimulateParticles<T, PARTICLETYPE> _sim;
+
+  std::deque<PARTICLETYPE<T> > _particles;
+
+  std::deque<PARTICLETYPE<T> > _shadowParticles;
+  std::list<std::shared_ptr<Force3D<T, PARTICLETYPE> > > _forces;
+  std::list<std::shared_ptr<Boundary3D<T, PARTICLETYPE> > > _boundaries;
+
+  std::vector<T> _physPos;
+  std::vector<T> _physExtend;
+
+
+  /// Integration methods, each need a special template particle
+  void velocityVerlet1(T dT);
+  void velocityVerlet2(T dT);
+//  void implicitEuler(T dT, AnalyticalF3D<T,T>& getvel);
+//  void adamBashforth4(T dT);
+//  void predictorCorrector1(T dT);
+//  void predictorCorrector2(T dT);
+  void rungeKutta4_1(T dt);
+  void rungeKutta4_2(T dt);
+  void rungeKutta4_3(T dt);
+  void rungeKutta4_4(T dt);
+  void rungeKutta4(T dT);
+  void updateParticleDistribution();
+};
+
+// Magnetic particle type
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::integrateTorqueMag(double dT);
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::integrateTorqueMag(double dT, std::set<int> sActivityOfParticle);
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::computeForce();
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::computeForce(std::set<int> sActivityOfParticle);
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::resetMag();
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::resetMag(std::set<int> sActivityOfParticle);
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::explicitEuler(double dT, bool scale);
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::explicitEuler(double dT, std::set<int> sActivityOfParticle, bool scale);
+
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::setOverlapZero();
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::setOverlapZeroForCombinationWithMechContactForce();
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::partialElasticImpact(double restitutionCoeff);
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::partialElasticImpactV2(double restitutionCoeff);
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::partialElasticImpactForCombinationWithMechContactForce(double restitutionCoeff);
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::findAgglomerates();
+template<>
+void ParticleSystem3D<double, MagneticParticle3D>::initAggloParticles();
+
+}  //namespace olb
+#endif /* PARTICLE_SYSTEM_3D_H */