/* 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
*
*
* 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
#include
#include
#include
#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
class SuperLatticeInterpPhysVelocity3D;
template class PARTICLETYPE>
class Force3D;
template class PARTICLETYPE>
class Boundary3D;
template class PARTICLETYPE>
class ParticleSystem3D;
template class PARTICLETYPE>
class SuperParticleSystem3D;
template class PARTICLETYPE>
class SuperParticleSysVtuWriter;
template
class SuperParticleSysVtuWriterMag;
template class PARTICLETYPE>
class SimulateParticles {
public:
SimulateParticles(ParticleSystem3D* 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 sActivityOfParticle, bool scale = false) {simulate(dT, scale);}
inline void simulateWithTwoWayCoupling_Mathias ( T dT,
ForwardCouplingModel& forwardCoupling,
BackCouplingModel& 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& forwardCoupling,
BackCouplingModel& 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* _pSys;
};
template
class SimulateParticles {
public:
SimulateParticles(ParticleSystem3D* ps)
: _pSys(ps)
{
}
inline void simulate(T dT, bool scale = false)
{
_pSys->computeForce();
_pSys->explicitEuler(dT, scale);
_pSys->integrateTorque(dT);
}
private:
ParticleSystem3D* _pSys;
};
template
class SimulateParticles {
public:
SimulateParticles(ParticleSystem3D* 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& forwardCoupling,
BackCouplingModel& 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& forwardCoupling,
BackCouplingModel& 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 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* _pSys;
};
template class PARTICLETYPE>
class ParticleSystem3D {
public:
/// Default constructor for ParticleSystem
ParticleSystem3D() = default;
/// Constructor for ParticleSystem
ParticleSystem3D(int iGeometry, SuperGeometry3D& superGeometry);
/// Copy constructor for ParticleSystem
ParticleSystem3D(const ParticleSystem3D& pS);
/// Move constructor for ParticleSystem
ParticleSystem3D(ParticleSystem3D && pS);
/// Destructor for ParticleSystem
virtual ~ParticleSystem3D() {}
virtual void simulate(T deltatime, bool scale = false);
virtual void simulateWithTwoWayCoupling_Mathias ( T dT,
ForwardCouplingModel& forwardCoupling,
BackCouplingModel& backCoupling,
int material, int subSteps, bool scale );
virtual void simulateWithTwoWayCoupling_Davide ( T dT,
ForwardCouplingModel& forwardCoupling,
BackCouplingModel& backCoupling,
int material, int subSteps, bool scale );
// multiple collision models
virtual void simulate(T deltatime, std::set 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& p);
/// Removes all particles from system
void clearParticles();
/// Add a force to ParticleSystem
void addForce(std::shared_ptr > pF);
/// Add a boundary to ParticleSystem
void addBoundary(std::shared_ptr > pF);
/// Get reference to a particle in the ParticleSystem
/// runs over all particles incl. shadow particles
PARTICLETYPE& operator[](const int i);
const PARTICLETYPE& operator[](const int i) const;
/// Set velocity of all particles to fluid velocity
template
void setVelToFluidVel(SuperLatticeInterpPhysVelocity3D &);
/// Set particle velocity to analytical velocity (e.g. as inital condition
void setVelToAnalyticalVel(AnalyticalConst3D&);
/// Set global coordinates and extends of Particlesystem (SI units)
void setPosExt(std::vector physPos, std::vector physExtend);
/// Get global coordinates and extends of Particlesystem (SI units)
const std::vector& getPhysPos();
const std::vector& 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 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 i, std::pair j) { return (i.second < j.second);}
} getMinDistPartObj;
void getMinDistParticle (std::vector> ret_matches);
/// Compute boundary contact
void computeBoundary();
/// Set boundary detection algorithm (for future features)
void setContactDetection(ContactDetection& contactDetection);
ContactDetection* getContactDetection();
/// Particle-Fluid interaction for subgrid scale particles
// template class DESCRIPTOR>
// void particleOnFluid(BlockLatticeStructure3D& bLattice,
// Cuboid3D& cuboid, int overlap, T eps,
// BlockGeometryStructure3D& bGeometry);
// template class DESCRIPTOR>
// void resetFluid(BlockLatticeStructure3D& bLattice,
// Cuboid3D& cuboid, int overlap);
friend class SuperParticleSystem3D;
friend class SuperParticleSysVtuWriter;
friend class SuperParticleSysVtuWriterMag;
friend class SimulateParticles;
//std::map 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 sActivityOfParticle, bool scale = false) {explicitEuler(dT, scale);};
bool executeForwardCoupling(ForwardCouplingModel& forwardCoupling);
bool executeBackwardCoupling(BackCouplingModel& backwardCoupling, int material, int subSteps=1);
ContactDetection* getDetection()
{
return _contactDetection;
}
int getIGeometry()
{
return _iGeometry;
}
/// returns deque of pointer to particles (not shadow particles)
/// contained in a particleSystem3D
std::deque*> getParticlesPointer();
/// returns deque of pointer to all particles (incl. shadow particles)
/// contained in a particleSystem3D
std::deque*> getAllParticlesPointer();
/// returns deque of pointer to all shadow particles
/// contained in a particleSystem3D
std::deque*> getShadowParticlesPointer();
/// returns deque of particles (no shadow particles)
/// contained in a particleSystem3D
std::deque> getParticles()
{
return _particles;
}
/// returns shared pointer of forces
std::list > > getForcesPointer();
/// Deque of Lists of agglomerated particles
std::deque*>> _Agglomerates;
protected:
void integrateTorque(T dT);
void integrateTorqueMag(T dT) {};
// multiple collision models
void integrateTorqueMag(T dT, std::set sActivityOfParticle) {};
void resetMag() {};
// multiple collision models
void resetMag(std::set 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& p);
mutable OstreamManager clout;
int _iGeometry = -1;
SuperGeometry3D& _superGeometry;
ContactDetection* _contactDetection;
SimulateParticles _sim;
std::deque > _particles;
std::deque > _shadowParticles;
std::list > > _forces;
std::list > > _boundaries;
std::vector _physPos;
std::vector _physExtend;
/// Integration methods, each need a special template particle
void velocityVerlet1(T dT);
void velocityVerlet2(T dT);
// void implicitEuler(T dT, AnalyticalF3D& 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::integrateTorqueMag(double dT);
template<>
void ParticleSystem3D::integrateTorqueMag(double dT, std::set sActivityOfParticle);
template<>
void ParticleSystem3D::computeForce();
template<>
void ParticleSystem3D::computeForce(std::set sActivityOfParticle);
template<>
void ParticleSystem3D::resetMag();
template<>
void ParticleSystem3D::resetMag(std::set sActivityOfParticle);
template<>
void ParticleSystem3D::explicitEuler(double dT, bool scale);
template<>
void ParticleSystem3D::explicitEuler(double dT, std::set sActivityOfParticle, bool scale);
template<>
void ParticleSystem3D::setOverlapZero();
template<>
void ParticleSystem3D::setOverlapZeroForCombinationWithMechContactForce();
template<>
void ParticleSystem3D::partialElasticImpact(double restitutionCoeff);
template<>
void ParticleSystem3D::partialElasticImpactV2(double restitutionCoeff);
template<>
void ParticleSystem3D::partialElasticImpactForCombinationWithMechContactForce(double restitutionCoeff);
template<>
void ParticleSystem3D::findAgglomerates();
template<>
void ParticleSystem3D::initAggloParticles();
} //namespace olb
#endif /* PARTICLE_SYSTEM_3D_H */