|
HG-MD
1
|
|
HG-MD
1
|
A class thast defines and solves a MD problem. More...
#include <MD.h>
Public Member Functions | |
| void | constructor () |
| A public constructor, which sets defaults so the problem can be solve off the shelf. | |
| MD () | |
| MD (STD_save &other) | |
| virtual | ~MD () |
| void | info () |
| Set up a virtual info this will be provided from the inhertiance. | |
| void | solve () |
| The work horse of the code. | |
| void | solve (unsigned int argc, char *argv[]) |
| Read arguments before solving. | |
| void | solveWithMDCLR () |
| Tries to solve the problem using MDCLR. | |
| Mdouble | get_t () |
| Access function for the time. | |
| void | set_t (Mdouble new_t) |
| Access function for the time. | |
| void | set_NWall (int N) |
| Allows the number of walls to be changed. | |
| int | get_NWall () |
| Allows the number of walls to be accessed. | |
| int | get_NSpecies () |
| Allows the number of Species to be accessed. | |
| vector< CWall > & | get_Walls (void) |
| Allows the walls to be copied. | |
| CWall & | get_Walls (int i) |
| Allows the walls to be accessed. | |
| vector< CWallPeriodic > & | get_WallsPeriodic (void) |
| Allows the periodic walls to be accessed. | |
| CWallPeriodic & | get_WallsPeriodic (int i) |
| Allows the periodic walls to be accessed. | |
| vector< CSpecies > & | get_Species (void) |
| Allows the species to be copied. | |
| CSpecies * | get_Species (int i) |
| Allows the species to be accessed. | |
| CSpecies * | get_MixedSpecies (int i, int j) |
| Allows the mixed species to be accessed. | |
| void | set_MixedSpecies (int i, int j, CSpecies &S) |
| Allows the mixed species to be set. | |
| void | set_NWallPeriodic (int N) |
| Allows the number of periodic walls to be changed. | |
| int | get_NWallPeriodic () |
| Allows the number of periodic walls to be accessed. | |
| void | set_tmax (Mdouble new_tmax) |
| Allows the upper time limit to be changed. | |
| Mdouble | get_tmax () |
| Allows the upper time limit to be accessed. | |
| ParticleHandler & | get_ParticleHandler () |
| void | set_savecount (int new_) |
| Allows the number of time steps between saves to be changed, see also set_number_of_saves. | |
| void | set_save_count_all (int new_) |
| void | set_save_count_data (int new_) |
| void | set_save_count_ene (int new_) |
| void | set_save_count_stat (int new_) |
| void | set_save_count_fstat (int new_) |
| int | get_savecount () |
| Allows the number of time steps between saves to be accessed. | |
| int | get_save_count () |
| int | get_save_count_data () |
| int | get_save_count_ene () |
| int | get_save_count_stat () |
| int | get_save_count_fstat () |
| void | set_do_stat_always (bool new_) |
| Sets how often the data is saved using the number of saves wanted, tmax, and dt. See also set_savecount. | |
| void | set_number_of_saves (Mdouble N) |
| void | set_number_of_saves_all (Mdouble N) |
| void | set_number_of_saves_data (Mdouble N) |
| void | set_number_of_saves_ene (Mdouble N) |
| void | set_number_of_saves_stat (Mdouble N) |
| void | set_number_of_saves_fstat (Mdouble N) |
| void | set_plastic_k1_k2max_kc_depth (Mdouble k1_, Mdouble k2max_, Mdouble kc_, Mdouble depth_, unsigned int indSpecies=0) |
| Allows the plastic constants to be changed. | |
| void | set_k1 (Mdouble new_, unsigned int indSpecies=0) |
| void | set_k2max (Mdouble new_, unsigned int indSpecies=0) |
| void | set_kc (Mdouble new_, unsigned int indSpecies=0) |
| void | set_depth (Mdouble new_, unsigned int indSpecies=0) |
| Mdouble | get_k1 (unsigned int indSpecies=0) |
| Allows the plastic constants to be accessed. | |
| Mdouble | get_k2max (unsigned int indSpecies=0) |
| Mdouble | get_kc (unsigned int indSpecies=0) |
| Mdouble | get_depth (unsigned int indSpecies=0) |
| Mdouble | get_plastic_dt (Mdouble mass, unsigned int indSpecies=0) |
| void | set_k (Mdouble new_, unsigned int indSpecies=0) |
| Allows the spring constant to be changed. | |
| Mdouble | get_k (int indSpecies=0) |
| Allows the spring constant to be accessed. | |
| void | set_kt (Mdouble new_, unsigned int indSpecies=0) |
| Allows the spring constant to be changed. | |
| Mdouble | get_kt (int indSpecies=0) |
| Allows the spring constant to be accessed. | |
| void | set_krolling (Mdouble new_, unsigned int indSpecies=0) |
| Allows the spring constant to be changed. | |
| Mdouble | get_krolling (int indSpecies=0) |
| Allows the spring constant to be accessed. | |
| void | set_ktorsion (Mdouble new_, unsigned int indSpecies=0) |
| Allows the spring constant to be changed. | |
| Mdouble | get_ktorsion (int indSpecies=0) |
| Allows the spring constant to be accessed. | |
| void | set_rho (Mdouble new_, unsigned int indSpecies=0) |
| Allows the density to be changed. | |
| Mdouble | get_rho (int indSpecies=0) |
| Allows the density to be accessed. | |
| void | set_dispt (Mdouble new_, unsigned int indSpecies=0) |
| Allows the tangential viscosity to be changed. | |
| Mdouble | get_dispt (unsigned int indSpecies=0) |
| Allows the tangential viscosity to be accessed. | |
| void | set_disprolling (Mdouble new_, unsigned int indSpecies=0) |
| Allows the tangential viscosity to be changed. | |
| Mdouble | get_disprolling (unsigned int indSpecies=0) |
| Allows the tangential viscosity to be accessed. | |
| void | set_disptorsion (Mdouble new_, unsigned int indSpecies=0) |
| Allows the tangential viscosity to be changed. | |
| Mdouble | get_disptorsion (unsigned int indSpecies=0) |
| Allows the tangential viscosity to be accessed. | |
| void | set_disp (Mdouble new_, unsigned int indSpecies=0) |
| Allows the normal dissipation to be changed. | |
| Mdouble | get_disp (unsigned int indSpecies=0) |
| Allows the normal dissipation to be accessed. | |
| void | set_dissipation (Mdouble new_, unsigned int indSpecies=0) |
| Allows the normal dissipation to be changed. | |
| Mdouble | get_dissipation (unsigned int indSpecies=0) |
| Allows the normal dissipation to be accessed. | |
| void | set_mu (Mdouble new_, unsigned int indSpecies=0) |
| Allows the Coulomb friction coefficient to be changed. | |
| Mdouble | get_mu (unsigned int indSpecies=0) |
| Allows the Coulomb friction coefficient to be accessed. | |
| void | set_murolling (Mdouble new_, unsigned int indSpecies=0) |
| Allows the Coulomb friction coefficient to be changed. | |
| Mdouble | get_murolling (unsigned int indSpecies=0) |
| Allows the Coulomb friction coefficient to be accessed. | |
| void | set_mutorsion (Mdouble new_, unsigned int indSpecies=0) |
| Allows the Coulomb friction coefficient to be changed. | |
| Mdouble | get_mutorsion (unsigned int indSpecies=0) |
| Allows the Coulomb friction coefficient to be accessed. | |
| void | set_dim_particle (int new_, unsigned int indSpecies=0) |
| Allows the dimension of the particle (f.e. for mass) to be changed. | |
| int | get_dim_particle (unsigned int indSpecies=0) |
| Allows the dimension of the particle (f.e. for mass) to be accessed. | |
| int | get_save_data_data () |
| Returns the data counter. | |
| int | get_save_data_ene () |
| int | get_save_data_fstat () |
| int | get_save_data_stat () |
| int | get_do_stat_always () |
| void | set_k_and_restitution_coefficient (Mdouble k_, Mdouble eps, Mdouble mass, unsigned int indSpecies=0) |
| Sets k, disp such that it matches a given tc and eps for a collision of two copies of P. | |
| void | set_collision_time_and_restitution_coefficient (Mdouble tc, Mdouble eps, Mdouble mass, unsigned int indSpecies=0) |
| Sets k, disp such that it matches a given tc and eps for a collision of two copies of P. | |
| void | set_collision_time_and_restitution_coefficient (Mdouble tc, Mdouble eps, Mdouble mass1, Mdouble mass2, unsigned int indSpecies=0) |
| Set k, disp such that is matches a given tc and eps for a collision of two different masses. | |
| void | set_collision_time_and_normal_and_tangential_restitution_coefficient (Mdouble tc, Mdouble eps, Mdouble beta, Mdouble mass1, Mdouble mass2, unsigned int indSpecies=0) |
| See CSpecies::set_collision_time_and_normal_and_tangential_restitution_coefficient. | |
| void | set_collision_time_and_normal_and_tangential_restitution_coefficient_nodispt (Mdouble tc, Mdouble eps, Mdouble beta, Mdouble mass1, Mdouble mass2, unsigned int indSpecies=0) |
| See CSpecies::set_collision_time_and_normal_and_tangential_restitution_coefficient. | |
| Mdouble | get_collision_time (Mdouble mass, unsigned int indSpecies=0) |
| Calculates collision time for two copies of a particle of given disp, k, mass. | |
| Mdouble | get_restitution_coefficient (Mdouble mass, unsigned int indSpecies=0) |
| Calculates restitution coefficient for two copies of given disp, k, mass. | |
| Mdouble | get_xmin () |
| Get xmin. | |
| Mdouble | get_xmax () |
| Get xmax. | |
| Mdouble | get_ymin () |
| Gets ymin. | |
| Mdouble | get_ymax () |
| Gets ymax. | |
| Mdouble | get_zmin () |
| Gets zmin. | |
| Mdouble | get_zmax () |
| Gets zmax. | |
| void | set_xmin (Mdouble new_xmin) |
| Sets xmin and walls, assuming the standard definition of walls as in the default constructor. | |
| void | set_ymin (Mdouble new_ymin) |
| void | set_zmin (Mdouble new_zmin) |
| Sets ymin and walls, assuming the standard definition of walls as in the default constructor. | |
| void | set_xmax (Mdouble new_xmax) |
| Sets xmax and walls, assuming the standard definition of walls as in the default constructor. | |
| void | set_ymax (Mdouble new_ymax) |
| Sets ymax and walls, assuming the standard definition of walls as in the default constructor. | |
| void | set_zmax (Mdouble new_zmax) |
| Sets ymax and walls, assuming the standard definition of walls as in the default constructor. | |
| void | set_dt (Mdouble new_dt) |
| Allows the time step dt to be changed. | |
| Mdouble | get_dt () |
| Allows the time step dt to be accessed. | |
| void | set_name (const char *name) |
| Sets the name of the problem, used for the same data files. | |
| void | set_xballs_colour_mode (int new_cmode) |
| Set the xball output mode. | |
| void | set_xballs_cmode (int new_cmode) |
| int | get_xballs_cmode () |
| void | set_xballs_vector_scale (double new_vscale) |
| Set the scale of vectors in xballs. | |
| double | get_xballs_vscale () |
| void | set_xballs_additional_arguments (string new_) |
| Set the additional arguments for xballs. | |
| string | get_xballs_additional_arguments () |
| void | set_xballs_scale (Mdouble new_scale) |
| Set the scale of the xballs problem. The default is fit to screen. | |
| double | get_xballs_scale () |
| void | set_gravity (Vec3D new_gravity) |
| Allows the gravitational acceleration to be changed. | |
| Vec3D | get_gravity () |
| Allows the gravitational acceleration to be accessed. | |
| void | set_dim (int new_dim) |
| Allows the dimension of the simulation to be changed. | |
| int | get_dim () |
| Allows the dimension of the simulation to be accessed. | |
| int | get_restart_version () |
| Gets restart_version. | |
| void | set_restart_version (int new_) |
| Sets restart_version. | |
| bool | get_restarted () |
| Gets restarted. | |
| Mdouble | get_max_radius () |
| Sets restarted. | |
| void | set_restarted (bool new_) |
| bool | get_append () |
| Gets restarted. | |
| void | set_append (bool new_) |
| Sets restarted. | |
| Mdouble | get_ene_ela () |
| Gets ene_ela. | |
| void | set_ene_ela (Mdouble new_) |
| Sets ene_ela. | |
| void | add_ene_ela (Mdouble new_) |
| Sets ene_ela. | |
| bool | get_Hertzian () |
| void | set_Hertzian (bool new_) |
| void | Remove_Particle (int IP) |
| This function removes partice IP from the vector of particles by moving the last particle in the vector to the position if IP Also it checks if the moved Particle has any tangentialsspring-information, which needs to be moved to a different particle, because tangential spring information always needs to be stored in the real particle with highest particle index. | |
| Mdouble | get_Mass_from_Radius (Mdouble radius, int indSpecies=0) |
| Mdouble | get_maximum_velocity (Particle &P) |
| Calculates the maximum velocity allowed for a collision of two copies of P (for higher velocities particles could pass through each other) | |
| virtual Particle * | get_SmallestParticle () |
| virtual Particle * | get_LargestParticle () |
| virtual void | removeParticle (int iP) |
| Mdouble | get_maximum_velocity () |
| void | set_dt_by_mass (Mdouble mass) |
| Sets dt to 1/50-th of the collision time for two particles of mass P. | |
| void | set_dt (Particle &P) |
| Sets dt to 1/50-th of the collision time for two copies of P. | |
| void | set_dt () |
| Sets dt to 1/50-th of the smallest possible collision time. | |
| virtual void | setup_particles_initial_conditions () |
| This function allows the initial conditions of the particles to be set, by default locations is random. | |
| virtual void | create_xballs_script () |
| This creates a scipt which can be used to load the xballs problem to display the data just generated. | |
| virtual double | getInfo (Particle &P UNUSED) |
| Allows the user to set what is written into the info column in the data file. | |
| virtual void | save_restart_data () |
| Stores all MD data. | |
| int | load_restart_data () |
| Loads all MD data. | |
| int | load_restart_data (string filename) |
| void | statistics_from_restart_data (const char *name) |
| Loads all MD data and plots statistics for all timesteps in the .data file. | |
| virtual void | write (std::ostream &os) |
| Writes all MD data. | |
| virtual void | read (std::istream &is) |
| Reads all MD data. | |
| virtual void | write_v1 (std::ostream &os) |
| Writes all MD data. | |
| virtual void | read_v1 (std::istream &is) |
| Reads all MD data. | |
| virtual void | read_v2 (std::istream &is) |
| bool | load_from_data_file (const char *filename, unsigned int format=0) |
| This allows particle data to be reloaded from data files. | |
| bool | load_par_ini_file (const char *filename) |
| allows the user to read par.ini files (useful to read MDCLR files) | |
| bool | read_next_from_data_file (unsigned int format=0) |
| int | read_dim_from_data_file () |
| bool | find_next_data_file (Mdouble tmin, bool verbose=true) |
| virtual void | print (std::ostream &os, bool print_all=false) |
| Outputs MD. | |
| void | add_Species (CSpecies &S) |
| void | add_Species (void) |
| void | set_format (int new_) |
| int | readArguments (unsigned int argc, char *argv[]) |
| Can interpret main function input arguments that are passed by the driver codes. | |
| virtual int | readNextArgument (unsigned int &i, unsigned int argc, char *argv[]) |
Public Attributes | |
| RNG | random |
Protected Member Functions | |
| virtual void | compute_all_forces () |
| This does the force computation. | |
| virtual void | compute_internal_forces (Particle *i) |
| Computes the forces between particles (internal in the sence that the sum over all these forces is zero i.e. fully modelled forces) | |
| virtual void | compute_internal_forces (Particle *P1, Particle *P2) |
| Computes the forces between particles (internal in the sence that the sum over all these forces is zero i.e. fully modelled forces) | |
| void | compute_plastic_internal_forces (Particle *P1, Particle *P2) |
| Computes plastic forces between particles. | |
| virtual void | compute_external_forces (Particle *PI) |
| This is were the computation of external forces takes place (e.g. gravity) | |
| virtual void | compute_walls (Particle *PI) |
| This is were the walls are. | |
| virtual void | actions_before_time_loop () |
| This is actions before the start of the main time loop. | |
| virtual void | HGRID_actions_before_time_loop () |
| This is actions before the start of the main time loop. | |
| virtual void | HGRID_actions_before_time_step () |
| This is action before the time step is started. | |
| virtual void | HGRID_InsertParticleToHgrid (Particle *obj UNUSED) |
| This is action before the time step is started. | |
| virtual void | HGRID_UpdateParticleInHgrid (Particle *obj UNUSED) |
| virtual void | HGRID_RemoveParticleFromHgrid (Particle *obj UNUSED) |
| virtual bool | get_HGRID_UpdateEachTimeStep () |
| virtual void | actions_before_time_step () |
| This is action before the time step is started. | |
| virtual void | actions_after_solve () |
| This is actions at the end of the code, but before the files are closed. | |
| virtual void | actions_after_time_step () |
| This is action after the time step is finished. | |
| virtual void | output_xballs_data () |
| Output xball data for Particle i. | |
| virtual void | output_xballs_data_particle (int i) |
| virtual void | start_ene () |
| Functions for ene file. | |
| virtual void | fstat_header () |
| virtual void | output_ene () |
| This function outputs statistical data - The default is to compute the rotational kinetic engergy, linear kinetic energy, and the centre of mass. | |
| virtual void | initialize_statistics () |
| Functions for statistics. | |
| virtual void | output_statistics () |
| virtual void | gather_statistics_collision (int index1 UNUSED, int index2 UNUSED, Vec3D Contact UNUSED, Mdouble delta UNUSED, Mdouble ctheta UNUSED, Mdouble fdotn UNUSED, Mdouble fdott UNUSED, Vec3D P1_P2_normal_ UNUSED, Vec3D P1_P2_tangential UNUSED) |
| virtual void | process_statistics (bool usethese UNUSED) |
| virtual void | finish_statistics () |
| virtual void | set_initial_pressures_for_pressure_controlled_walls () |
| virtual void | do_integration_before_force_computation (Particle *pI) |
| This is were the integration is done. | |
| virtual void | HGRID_update_move (Particle *, Mdouble) |
| virtual void | HGRID_actions_before_integration () |
| virtual void | HGRID_actions_after_integration () |
| virtual void | do_integration_after_force_computation (Particle *pI) |
| This is were the integration is done. | |
| virtual void | InitBroadPhase () |
| Initialisation of Broad Phase Information (Default no Broad Phase so empty) | |
| virtual void | broad_phase (Particle *i) |
| Broad phase of contact detection goes here. Default check all contacts. | |
| void | set_FixedParticles (unsigned int n) |
| void | initialize_tangential_springs () |
| void | compute_particle_masses () |
| Computes the mass of each particle. | |
| virtual void | cout_time () |
| Couts time. | |
| virtual bool | continue_solve () |
| void | reset_DeltaMax () |
| sets the history parameter DeltaMax of all particles to zero | |
| void | reset_TangentialSprings () |
| sets the history parameter TangentialSprings of all particles to zero | |
Protected Attributes | |
| vector< CWall > | Walls |
| This is a particle class, which stores all information about each particle, i.e. | |
| vector< CWallPeriodic > | WallsPeriodic |
| vector< CSpecies > | Species |
| These are the particle parameters like dissipation etc. | |
Private Member Functions | |
| void | Remove_Duplicate_Periodic_Particles () |
| Remove periodic duplicate Particles (i.e. removes particles created by Check_and_Duplicate_Periodic_Particle(int i, int nWallPeriodic)). Note that between these two functions it is not allowed to create additional functions. | |
| virtual int | Check_and_Duplicate_Periodic_Particle (Particle *i, int nWallPeriodic) |
| Checks if Particle with index i (i.e. Particle[i]) is close to a periodic wall with index NWallPeriodic or higher. If so it creates a periodic duplicate of the particle at the end of the particle vector wich has an variable periodicFromParticle which is the index of the originating (real) particle. Then it calls itself recusrively to also get doubly Periodic shifted images. It returns the total number of Periodic Images created. | |
| int | Check_and_Duplicate_Periodic_Particles () |
| Calls Check_and_Duplicate_Periodic_Particle for all Particles curently in Particles[] and returns the number of particles created. | |
Private Attributes | |
| int | dim |
| The dimension of the simulation. | |
| Vec3D | gravity |
| Gravitational acceleration. | |
| Mdouble | max_radius |
| Mdouble | xmin |
| These store the size of the domain, assume walls at the ends. | |
| Mdouble | xmax |
| Mdouble | ymin |
| Mdouble | ymax |
| Mdouble | zmin |
| Mdouble | zmax |
| Mdouble | dt |
| These are numerical constants like the time step size. | |
| Mdouble | tmax |
| int | save_count_data |
| These are informations for saving. | |
| int | save_count_ene |
| int | save_count_stat |
| int | save_count_fstat |
| bool | save_data_data |
| bool | save_data_ene |
| bool | save_data_fstat |
| bool | save_data_stat |
| bool | do_stat_always |
| Mdouble | t |
| This stores the current time. | |
| Mdouble | ene_ela |
| used in force calculations | |
| int | xballs_cmode |
| Mdouble | xballs_vscale |
| Mdouble | xballs_scale |
| string | xballs_additional_arguments |
| int | format |
| int | restart_version |
| bool | restarted |
| bool | Hertzian |
| int | data_FixedParticles |
| ParticleHandler | particleHandler |
| bool | append |
| bool | rotation |
| int | PeriodicCreated |
| int | save_restart_data_counter |
Friends | |
| std::ostream & | operator<< (std::ostream &os, MD &md) |
A class thast defines and solves a MD problem.
| MD::MD | ( | ) | [inline] |
{
constructor();
#ifdef CONSTUCTOR_OUTPUT
cerr << "MD() finished"<<endl;
#endif
}
: STD_save(other) { constructor(); #ifdef CONSTUCTOR_OUTPUT cerr << "MD(STD_save& other) finished " << endl; #endif }
| virtual void MD::actions_after_solve | ( | ) | [inline, protected, virtual] |
| virtual void MD::actions_after_time_step | ( | ) | [inline, protected, virtual] |
This is action after the time step is finished.
Referenced by solve(), and statistics_from_restart_data().
{};
| virtual void MD::actions_before_time_loop | ( | ) | [inline, protected, virtual] |
This is actions before the start of the main time loop.
todo thomas: can this be erased?
Referenced by solve(), and statistics_from_restart_data().
| virtual void MD::actions_before_time_step | ( | ) | [inline, protected, virtual] |
This is action before the time step is started.
Reimplemented in Chute, and ChuteBottom.
Referenced by solve(), and statistics_from_restart_data().
{};
| void MD::add_ene_ela | ( | Mdouble | new_ | ) | [inline] |
Sets ene_ela.
{ene_ela+=new_;}
| void MD::add_Species | ( | CSpecies & | S | ) |
| void MD::add_Species | ( | void | ) | [inline] |
| virtual void MD::broad_phase | ( | Particle * | i | ) | [inline, protected, virtual] |
Broad phase of contact detection goes here. Default check all contacts.
Reimplemented in HGRID_base.
Referenced by compute_internal_forces().
{
for (vector<Particle*>::iterator it = particleHandler.begin(); (*it)!=i; ++it)
{
compute_internal_forces(i,*it);
}
}
| int MD::Check_and_Duplicate_Periodic_Particle | ( | Particle * | i, |
| int | nWallPeriodic | ||
| ) | [private, virtual] |
Checks if Particle with index i (i.e. Particle[i]) is close to a periodic wall with index NWallPeriodic or higher. If so it creates a periodic duplicate of the particle at the end of the particle vector wich has an variable periodicFromParticle which is the index of the originating (real) particle. Then it calls itself recusrively to also get doubly Periodic shifted images. It returns the total number of Periodic Images created.
References ParticleHandler::back(), ParticleHandler::copyAndAddParticle(), get_NWallPeriodic(), Particle::get_PeriodicFromParticle(), Particle::get_Position(), Particle::get_Radius(), HGRID_InsertParticleToHgrid(), HGRID_UpdateParticleInHgrid(), max_radius, particleHandler, Particle::set_periodicFromParticle(), and WallsPeriodic.
Referenced by Check_and_Duplicate_Periodic_Particles().
{
int C=0; //Number of particles created
for (int k=nWallPeriodic; k<get_NWallPeriodic(); k++) { //Loop over all still posible walls
if (WallsPeriodic[k].distance(*i)<i->get_Radius()+max_radius)
{
Particle F0=*i;
WallsPeriodic[k].shift_position(F0.get_Position());
//If Particle is Mdouble shifted, get correct original particle
Particle* From=i;
while(From->get_PeriodicFromParticle()!=NULL)
From=From->get_PeriodicFromParticle();
F0.set_periodicFromParticle(From);
particleHandler.copyAndAddParticle(F0);
HGRID_InsertParticleToHgrid(particleHandler.back());
HGRID_UpdateParticleInHgrid(particleHandler.back());
C++;
//Check for Mdouble shifted particles
C+=Check_and_Duplicate_Periodic_Particle(particleHandler.back(), k+1);
}
}
return(C);
}
| int MD::Check_and_Duplicate_Periodic_Particles | ( | ) | [private] |
Calls Check_and_Duplicate_Periodic_Particle for all Particles curently in Particles[] and returns the number of particles created.
References Check_and_Duplicate_Periodic_Particle(), ParticleHandler::get_NumberOfParticles(), ParticleHandler::get_Particle(), and particleHandler.
Referenced by solve(), and statistics_from_restart_data().
{
int N=particleHandler.get_NumberOfParticles(); //Required because number of particles increaes
int C=0;
for (int i=0; i<N; i++)
C+=Check_and_Duplicate_Periodic_Particle(particleHandler.get_Particle(i),0);
return(C);
}
| void MD::compute_all_forces | ( | ) | [protected, virtual] |
This does the force computation.
Reset all forces to zero
Now loop over all particles contacts computing force contributions
Now loop over all other particles looking for contacts
References ParticleHandler::begin(), compute_external_forces(), compute_internal_forces(), ParticleHandler::end(), and particleHandler.
Referenced by solve(), and statistics_from_restart_data().
{
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it)
{
(*it)->get_Force().set_zero();
(*it)->get_Torque().set_zero();
} //end reset forces loop
#ifdef DEBUG_OUTPUT
cerr << "Have all forces set to zero " << endl;
#endif
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it)
{
compute_internal_forces(*it);
//end inner loop over contacts.
compute_external_forces(*it);
}
//end outer loop over contacts.
}
| void MD::compute_external_forces | ( | Particle * | CI | ) | [protected, virtual] |
This is were the computation of external forces takes place (e.g. gravity)
This computes the external forces e.g.
here it is gravity and walls
Now add on gravity
Finally walls
References Particle::add_Force(), compute_walls(), Particle::get_Mass(), gravity, and Particle::isFixed().
Referenced by compute_all_forces().
{
CI->add_Force(gravity * CI->get_Mass());
if (!CI->isFixed()) compute_walls(CI);
}
| void MD::compute_internal_forces | ( | Particle * | i | ) | [protected, virtual] |
Computes the forces between particles (internal in the sence that the sum over all these forces is zero i.e. fully modelled forces)
References broad_phase().
Referenced by HGRID_2D::CheckCell(), HGRID_3D::CheckCell(), HGRID_2D::CheckCell_current(), HGRID_3D::CheckCell_current(), and compute_all_forces().
{
broad_phase(i);
}
| void MD::compute_internal_forces | ( | Particle * | P1, |
| Particle * | P2 | ||
| ) | [protected, virtual] |
Computes the forces between particles (internal in the sence that the sum over all these forces is zero i.e. fully modelled forces)
This computer the internal forces (internal in the sence that they sum to zero) i.e.
the fully modelled forces.
Tangential spring information is always store in the real particle with highest index When a Periodic contact is encountered it is always encoutered twice, but only applied when the real particle has the highest index of both real indices When a Particle is removed the tangential spring information has to be moved
Both options are up to first order the same (the first one is nicer becaus it always keeps the spring tangential, whereas the second one is in a nicer intergration form)
The second particle (i.e. the particle the force acts on) is always a flow particle
References Particle::add_Force(), Particle::add_Torque(), CTangentialSpring::delta, CSpecies::disp, CSpecies::disprolling, CSpecies::dispt, CSpecies::disptorsion, do_stat_always, dt, ene_ela, STD_save::fstat_file, gather_statistics_collision(), Particle::get_AngularVelocity(), get_Hertzian(), Particle::get_Index(), Particle::get_IndSpecies(), get_MixedSpecies(), Particle::get_PeriodicFromParticle(), Particle::get_Position(), Particle::get_Radius(), Particle::get_TangentialSprings(), Particle::get_Velocity(), Vec3D::GetLength(), Vec3D::GetLength2, Particle::isFixed(), CSpecies::k, CSpecies::krolling, CSpecies::kt, CSpecies::ktorsion, CSpecies::mu, CSpecies::murolling, CSpecies::mus, CSpecies::musrolling, CSpecies::mustorsion, CSpecies::mutorsion, CTangentialSpring::RollingSpring, save_data_ene, save_data_fstat, save_data_stat, CTangentialSprings::select_particle_spring(), Vec3D::set_zero(), CTangentialSpring::sliding, Species, sqr, t, and CTangentialSpring::TorsionSpring.
{
//this is because the rough bottom allows overlapping fixed particles
if (P1->isFixed()&&P2->isFixed()) return;
//Cases (start from P1 and P2 and go to PI and PJ
//1 Normal-Normal ->PI=max(P1,P2), PJ=min(P1,P2)
//2 Periodic-Normal ->if(P2>Real(P1)) (PI=P2 PJ=real(P1)) otherwise do nothing
//3 Normal-Periodic ->if(P1>Real(P2)) (PI=P1 PJ=real(P2)) otherwise do nothing
//4 Periodic-Periodic ->do nothing
//Just some statements to handle the 4 cases
Particle *PI, *PJ,*PJreal;
Particle *P1Per=P1->get_PeriodicFromParticle();
Particle *P2Per=P2->get_PeriodicFromParticle();
if(P1Per==NULL)
{
if(P2Per==NULL)
//N-N
if(P1->get_Index()<P2->get_Index())
{PI=P2;PJ=P1;PJreal=PJ;}
else
{PI=P1;PJ=P2;PJreal=PJ;}
else
//N-P
if(P1->get_Index()>P2Per->get_Index())
{PI=P1;PJ=P2;PJreal=P2Per;}
else return;
} else {
if(P2Per==NULL)
//P-N
if(P2->get_Index()>P1Per->get_Index())
{PI=P2;PJ=P1;PJreal=P1Per;}
else return;
else return;
}
Vec3D vrel, vrelt, forcet = Vec3D(0.0, 0.0, 0.0);
Vec3D deltat = Vec3D(0.0, 0.0, 0.0);
#ifdef DEBUG_OUTPUT
cerr << "In computing interal forces between particle "<<PI->get_Position()<<" and "<<PJ->get_Position()<<endl;
#endif
//Get the square of the distance between particle i and particle j
Mdouble dist_squared=GetDistance2(PI->get_Position(), PJ->get_Position());
Mdouble radii_sum=PI->get_Radius()+PJ->get_Radius();
#ifdef DEBUG_OUTPUT_FULL
cerr << "Square of distance between " << dist_squared << " square sum of radii " << radii_sum*radii_sum <<endl;
#endif
// True if the particles are in contact
if (dist_squared<(radii_sum*radii_sum))
{
// For particles of the same species, set species vector to Species(PI);
// for particles of different species, set species vector to MixedSpecies(PI,PJ)
CSpecies* pSpecies = (PI->get_IndSpecies()==PJ->get_IndSpecies())?&Species[PI->get_IndSpecies()]:get_MixedSpecies(PI->get_IndSpecies(),PJ->get_IndSpecies());
// Calculate distance between the particles
Mdouble dist=sqrt(dist_squared);
// Compute normal vector
Vec3D normal=(PI->get_Position()-PJ->get_Position())/dist;
// Compute the overlap between the particles
Mdouble deltan = radii_sum-dist;
// Compute the relative velocity vector v_ij
Vec3D vrel;
if (!pSpecies->mu) {
vrel=(PI->get_Velocity()-PJ->get_Velocity());
} else {
vrel=(PI->get_Velocity()-PJ->get_Velocity()) + Cross(normal, PI->get_AngularVelocity() * (PI->get_Radius() - .5 * deltan) + PJ->get_AngularVelocity() * (PJ->get_Radius() - .5 * deltan) );
}
// Compute the projection of vrel onto the normal (can be negative)
Mdouble vdotn=-Dot(vrel,normal);
// Compute normal force on particle i due to contact
Mdouble fdotn = pSpecies->k*deltan+pSpecies->disp*vdotn;
Vec3D force = normal * fdotn;
//If tangential forces are present
CTangentialSpring* TangentialSpring = NULL;
Vec3D forcet; forcet.set_zero();
Vec3D forcerolling; forcerolling.set_zero();
Vec3D forcetorsion; forcetorsion.set_zero();
if (pSpecies->mu || pSpecies->murolling || pSpecies->mutorsion) {
//call tangential spring
if (pSpecies->kt || pSpecies->krolling || pSpecies->ktorsion)
TangentialSpring = PI->get_TangentialSprings().select_particle_spring(PJreal->get_Index(), t, dt);
//Compute norm of normal force
Mdouble norm_fn = abs(fdotn);
//calculate sliding friction
if (pSpecies->mu) {
//Compute the tangential component of vrel
Vec3D vrelt=vrel+normal*vdotn;
//Compute norm of vrelt
Mdouble vdott=vrelt.GetLength();
if (pSpecies->kt) {
Vec3D* delta = &(TangentialSpring->delta);
//Integrate the spring
//(*delta) += vrelt * dt - Dot(*delta,normal)*normal;
(*delta) += (vrelt - Dot(*delta,PI->get_Velocity()-PJ->get_Velocity())*normal/dist) * dt;
//Calculate test force including viscous force
forcet = (-pSpecies->dispt) * vrelt - pSpecies->kt * (*delta);
Mdouble forcet2 = forcet.GetLength2();
//tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking),
//but the force is limited by Coulomb friction (sliding):
//f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
if( forcet2 <= sqr(pSpecies->mus*norm_fn) ) {
//sticking++;
TangentialSpring->sliding=false;
} else {
//sliding++;
TangentialSpring->sliding=true;
Mdouble norm_forcet = sqrt(forcet2);
forcet *= pSpecies->mu * norm_fn / norm_forcet;
(*delta) = (forcet + pSpecies->dispt * vrelt)/(-pSpecies->kt);
//~ forcet = (-pSpecies->mu * norm_fn) * GetUnitVector(*delta);
//~ (*delta) = forcet/(-pSpecies->kt);
}
//Add tangential force to total force
force += forcet;
} else { //if no tangential spring
//tangential forces are modelled by a damper of viscosity dispt (sticking),
//but the force is limited by Coulomb friction (sliding):
//f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
if (vdott*pSpecies->dispt <= pSpecies->mus*norm_fn) { //sticking++;
forcet = -pSpecies->dispt * vrelt;
} else { //sliding++;
//set force to Coulomb limit
forcet = -(pSpecies->mu * norm_fn / vdott) * vrelt;
}
//Add tangential force to total force
force += forcet;
}
}
//calculate rolling friction
if (pSpecies->murolling) {
//From Luding2008, objective rolling velocity (eq 15) w/o 2.0!
Mdouble reducedRadiusI = PI->get_Radius() - .5 * deltan;
Mdouble reducedRadiusJ = PJ->get_Radius() - .5 * deltan;
Mdouble reducedRadiusIJ = 2.0*reducedRadiusI*reducedRadiusJ/(reducedRadiusI+reducedRadiusJ);
Vec3D vrolling=-reducedRadiusIJ * Cross(normal, PI->get_AngularVelocity() - PJ->get_AngularVelocity());
if (pSpecies->krolling) {
Vec3D* RollingSpring = &(TangentialSpring->RollingSpring);
//Integrate the spring
(*RollingSpring) += vrolling * dt;// - Dot(*RollingSpring,normal)*normal;
//Calculate test force including viscous force
forcerolling = (-pSpecies->disprolling) * vrolling - pSpecies->krolling * (*RollingSpring);
Mdouble forcerolling2 = forcerolling.GetLength2();
//tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking),
//but the force is limited by Coulomb friction (sliding):
//f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
if( forcerolling2 <= sqr(pSpecies->musrolling*norm_fn) ) {
//sticking++;
} else {
//sliding++;
forcerolling *= pSpecies->murolling * norm_fn / sqrt(forcerolling2);
(*RollingSpring) = (forcerolling + pSpecies->disprolling * vrolling)/(-pSpecies->krolling);
}
//Add tangential force to torque
Vec3D Torque = reducedRadiusIJ * Cross(normal, forcerolling);
PI->add_Torque(Torque);
PJreal->add_Torque(-Torque);
}
}
//calculate torsive friction
if (pSpecies->mutorsion) {
//From Luding2008, spin velocity (eq 16) w/o 2.0!
Mdouble RadiusIJ = 2.0*PI->get_Radius()*PJ->get_Radius()/(PI->get_Radius()+PJ->get_Radius());
Vec3D vtorsion=RadiusIJ*Dot(normal,PI->get_AngularVelocity() - PJ->get_AngularVelocity())*normal;
if (pSpecies->ktorsion) {
//~ cout << "Error; not yet implemented" << endl;
//~ exit(-1);
Vec3D* TorsionSpring = &(TangentialSpring->TorsionSpring);
//Integrate the spring
(*TorsionSpring) = Dot((*TorsionSpring) + vtorsion * dt,normal)*normal;
//Calculate test force including viscous force
forcetorsion = (-pSpecies->disptorsion) * vtorsion - pSpecies->ktorsion * (*TorsionSpring);
Mdouble forcetorsion2 = forcetorsion.GetLength2();
//tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking),
//but the force is limited by Coulomb friction (sliding):
//f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
if( forcetorsion2 <= sqr(pSpecies->mustorsion*norm_fn) ) {
//sticking++;
} else {
//sliding++;
cout << "sliding " << endl;
forcetorsion *= pSpecies->mutorsion * norm_fn / sqrt(forcetorsion2);
(*TorsionSpring) = (forcetorsion + pSpecies->disptorsion * vtorsion)/(-pSpecies->ktorsion);
}
//Add tangential force to torque
Vec3D torque = RadiusIJ * forcetorsion;
PI->add_Torque(torque);
PJreal->add_Torque(-torque);
}
}
} //end if tangential forces
//Make force Hertzian (note: deltan is normalized by the normal distanct of two particles in contact, as in Silbert)
if (get_Hertzian()) force *= sqrt(deltan/(PI->get_Radius()+PJ->get_Radius()));
//Add forces to total force
PI->add_Force(force);
PJreal->add_Force(-force);
// Add torque due to tangential forces: t = Cross(l,f), l=dist*Wall.normal
if (pSpecies->mu) {
Vec3D cross = Cross(normal, force);
PI ->add_Torque(-cross * (PI->get_Radius() - .5 * deltan));
PJreal->add_Torque(-cross * (PJ->get_Radius() - .5 * deltan));
}
// output for ene and stat files:
if (save_data_ene) {
ene_ela += 0.5 * (pSpecies->k * sqr(deltan) + (TangentialSpring?
(pSpecies->kt * TangentialSpring->delta.GetLength2()
+pSpecies->krolling * TangentialSpring->RollingSpring.GetLength2()
+pSpecies->ktorsion * TangentialSpring->TorsionSpring.GetLength2()):0.0));
}
if (save_data_fstat||save_data_stat||do_stat_always) {
Mdouble fdott = forcet.GetLength();
Mdouble deltat_norm = TangentialSpring?(-TangentialSpring->delta.GetLength()):0.0;
Vec3D centre = 0.5 * (PI->get_Position() + PJ->get_Position());
if (!PI->isFixed())
{
if(save_data_stat||do_stat_always) gather_statistics_collision(PJreal->get_Index(),PI->get_Index(), centre, deltan, deltat_norm,fdotn,fdott,-normal,-(fdott?forcet/fdott:forcet));
if(save_data_fstat) fstat_file
<< t << " "
<< PJreal->get_Index() << " "
<< PI->get_Index() << " "
<< centre << " "
<< deltan << " "
<< deltat_norm << " "
<< fdotn << " "
<< fdott << " "
<< -normal << " "
<< -(fdott?forcet/fdott:forcet) << endl;
}
if (!PJreal->isFixed())
{
if(save_data_stat||do_stat_always) gather_statistics_collision(PI->get_Index(),PJreal->get_Index(), centre, deltan, deltat_norm,fdotn,fdott,normal,(fdott?forcet/fdott:forcet));
if(save_data_fstat) fstat_file
<< t << " "
<< PI->get_Index() << " "
<< PJreal->get_Index() << " "
<< centre << " "
<< deltan << " "
<< deltat_norm << " "
<< fdotn << " "
<< fdott << " "
<< normal << " "
<< (fdott?forcet/fdott:forcet) << endl;
}
}
} // end if particle i and j are overlapping
}
| void MD::compute_particle_masses | ( | ) | [inline, protected] |
Computes the mass of each particle.
Referenced by solve(), and statistics_from_restart_data().
{for (unsigned int i=0;i<particleHandler.get_NumberOfParticles();i++) particleHandler.get_Particle(i)->compute_particle_mass(Species);}
| void MD::compute_plastic_internal_forces | ( | Particle * | P1, |
| Particle * | P2 | ||
| ) | [protected] |
Computes plastic forces between particles.
This computer the internal forces for the plastic model.
Tangential spring information is always store in the real particle with highest index When a Periodic contact is encountered it is always encoutered twice, but only applied when the real particle has the highest index of both real indices When a Particle is removed the tangential spring information has to be moved
newcode begin
Both options are up to first order the same (the first one is nicer becaus it always keeps the spring tangential, whereas the second one is in a nicer intergration form)
The second particle (i.e. the particle the force acts on) is always a flow particle
References Particle::add_Force(), Particle::add_Torque(), CTangentialSpring::delta, CSpecies::depth, CSpecies::disp, CSpecies::dispt, do_stat_always, dt, ene_ela, STD_save::fstat_file, gather_statistics_collision(), Particle::get_AngularVelocity(), Particle::get_DeltaMaxs(), Particle::get_Index(), Particle::get_IndSpecies(), get_MixedSpecies(), Particle::get_PeriodicFromParticle(), Particle::get_Position(), Particle::get_Radius(), Particle::get_TangentialSprings(), Particle::get_Velocity(), Vec3D::GetLength(), Vec3D::GetLength2, Hertzian, Particle::isFixed(), CSpecies::k, CSpecies::k2max, CSpecies::kc, CSpecies::kt, CSpecies::mu, CSpecies::mus, save_data_ene, save_data_fstat, save_data_stat, CDeltaMaxs::select_particle(), CTangentialSprings::select_particle_spring(), Vec3D::set_zero(), Species, sqr, and t.
{
//this is because the rough bottom allows overlapping fixed particles
if (P1->isFixed()&&P2->isFixed()) return;
//Cases (start from P1 and P2 and go to PI and PJ
//1 Normal-Normal ->PI=max(P1,P2), PJ=min(P1,P2)
//2 Periodic-Normal ->if(P2>Real(P1)) (PI=P2 PJ=real(P1)) otherwise do nothing
//3 Normal-Periodic ->if(P1>Real(P2)) (PI=P1 PJ=real(P2)) otherwise do nothing
//4 Periodic-Periodic ->do nothing
//Just some statements to handle the 4 cases
Particle *PI, *PJ,*PJreal;
Particle *P1Per=P1->get_PeriodicFromParticle();
Particle *P2Per=P2->get_PeriodicFromParticle();
if(P1Per==NULL)
{
if(P2Per==NULL)
//N-N
if(P1->get_Index()<P2->get_Index())
{PI=P2;PJ=P1;PJreal=PJ;}
else
{PI=P1;PJ=P2;PJreal=PJ;}
else
//N-P
if(P1->get_Index()>P2Per->get_Index())
{PI=P1;PJ=P2;PJreal=P2Per;}
else return;
} else {
if(P2Per==NULL)
//P-N
if(P2->get_Index()>P1Per->get_Index())
{PI=P2;PJ=P1;PJreal=P1Per;}
else return;
else return;
}
#ifdef DEBUG_OUTPUT
cerr << "In computing interal forces between particle "<<PI->get_Position()<<" and "<<PJ->get_Position()<<endl;
#endif
//Get the square of the distance between particle i and particle j
Mdouble dist_squared=GetDistance2(PI->get_Position(), PJ->get_Position());
Mdouble radii_sum=PI->get_Radius()+PJ->get_Radius();
#ifdef DEBUG_OUTPUT_FULL
cerr << "Square of distance between " << dist_squared << " square sum of radii " << radii_sum*radii_sum <<endl;
#endif
// True if the particles are in contact
if (dist_squared<(radii_sum*radii_sum))
{
// For particles of the same species, set species vector to Species(PI);
// for particles of different species, set species vector to MixedSpecies(PI,PJ)
CSpecies* pSpecies = (PI->get_IndSpecies()==PJ->get_IndSpecies())?&Species[PI->get_IndSpecies()]:get_MixedSpecies(PI->get_IndSpecies(),PJ->get_IndSpecies());
// Calculate distance between the particles
Mdouble dist=sqrt(dist_squared);
// Compute normal vector
Vec3D normal=(PI->get_Position()-PJ->get_Position())/dist;
// Compute the overlap between the particles
Mdouble deltan = radii_sum-dist;
// Compute the relative velocity vector v_ij
Vec3D vrel;
if (!pSpecies->mu) {
vrel=(PI->get_Velocity()-PJ->get_Velocity());
} else {
vrel=(PI->get_Velocity()-PJ->get_Velocity()) + Cross(normal, PI->get_AngularVelocity() * (PI->get_Radius() - .5 * deltan) + PJ->get_AngularVelocity() * (PJ->get_Radius() - .5 * deltan) );
}
// BEGIN add viscous normal force
// Compute the projection of vrel onto the normal (can be negative)
Mdouble vdotn=-Dot(vrel,normal);
// Compute normal force on particle i due to contact
Mdouble fdotn = pSpecies->disp*vdotn;
// END add viscous normal force
// BEGIN add elastic force
//calculate deltastar, the overlap above which k2max becomes active (the 'fluid branch')
Mdouble deltastar = (pSpecies->k2max/(pSpecies->k2max-pSpecies->k))*pSpecies->depth*((2*PI->get_Radius()*PJ->get_Radius())/(PI->get_Radius()+PJ->get_Radius()));
//retrieve history parameter deltamax, the max. achieved overlap
Mdouble *deltamax = PI->get_DeltaMaxs().select_particle(PJreal->get_Index(), t, dt);
//update deltastar
*deltamax = min(deltastar,max(*deltamax,deltan));
//calculate k2(deltamax)
Mdouble k2;
if (deltan>deltastar) {
k2 = pSpecies->k2max;
} else {
k2 = pSpecies->k+(pSpecies->k2max-pSpecies->k)*((*deltamax)/deltastar);
}
//calculate delta0(deltamax), the overlap where the force is zero
Mdouble delta0 = (k2-pSpecies->k)/k2*(*deltamax);
//add elastic force
//cout << k2*(deltan-(delta0))-pSpecies->k*deltan << endl;
if (deltan>deltastar) {
fdotn += pSpecies->k2max*(deltan-(delta0));
} else if (k2*(deltan-(delta0))>=pSpecies->k*deltan){
fdotn += pSpecies->k*deltan;
} else if (k2*(deltan-delta0)>=-pSpecies->kc*deltan){
fdotn += k2*(deltan-delta0);
} else {
fdotn += -pSpecies->kc*deltan;
*deltamax=(k2+pSpecies->kc)/(k2-pSpecies->k)*deltan;
}
Vec3D force = normal * fdotn;
// END add elastic force
//If tangential forces are present
Vec3D forcet; forcet.set_zero();
Vec3D deltat; deltat.set_zero();
if (pSpecies->mu) {
//Compute the tangential component of vrel
Vec3D vrelt=vrel+normal*vdotn;
//Compute norm of vrelt
Mdouble vdott=vrelt.GetLength();
//Compute norm of normal force
Mdouble norm_fn = abs(fdotn);
if (!pSpecies->kt) { //if no tangential spring
//tangential forces are modelled by a damper of viscosity dispt (sticking),
//but the force is limited by Coulomb friction (sliding):
//f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
if (vdott*pSpecies->dispt <= pSpecies->mus*norm_fn) { //sticking++;
forcet = -pSpecies->dispt * vrelt;
} else { //sliding++;
//set force to Coulomb limit
forcet = -(pSpecies->mu * norm_fn / vdott) * vrelt;
}
} else { //if with tangential spring
//Retrieve the spring
CTangentialSpring* TangentialSpring = PI->get_TangentialSprings().select_particle_spring(PJreal->get_Index(), t, dt);
Vec3D* delta = &(TangentialSpring->delta);
//Integrate the spring
//(*delta) += vrelt * dt - Dot(*delta,normal)*normal;
(*delta) += (vrelt - Dot(*delta,PI->get_Velocity()-PJ->get_Velocity())*normal/dist) * dt;
deltat = (*delta);
//Calculate test force including viscous force
forcet = (-pSpecies->dispt) * vrelt - pSpecies->kt * deltat;
Mdouble forcet2 = forcet.GetLength2();
//tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking),
//but the force is limited by Coulomb friction (sliding):
//f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
if( forcet2 <= sqr(pSpecies->mus*norm_fn) ) {
//sticking++;
} else {
//sliding++;
Mdouble norm_forcet = sqrt(forcet2);
forcet *= pSpecies->mu * norm_fn / norm_forcet;
(*delta) = -(forcet + pSpecies->dispt * vrelt)/pSpecies->kt;
//This line should be removed before release it is the old tangential model (the new is shown to be better).
//(*delta) = forcet/(-pSpecies->kt);
}
} //end if tangential spring
//Add tangential force to total force
force += forcet;
} //end if tangential forces
//Make force Hertzian (note: deltan is normalized by the normal distanct of two particles in contact, as in Silbert)
if (Hertzian) force *= sqrt(deltan/(PI->get_Radius()+PJ->get_Radius()));
PI ->add_Force( force);
PJreal->add_Force(-force);
// Add torque due to tangential forces: t = Cross(l,f), l=dist*Wall.normal
if (pSpecies->mu) {
Vec3D cross = -Cross(normal, force);
PI ->add_Torque(cross * (PI->get_Radius() - .5 * deltan));
PJreal->add_Torque(cross * (PI->get_Radius() - .5 * deltan));
}
// output for ene and stat files:
if (save_data_ene) {
ene_ela += 0.5 * (pSpecies->k * sqr(deltan) + pSpecies->kt * deltat.GetLength2());
}
if (save_data_fstat||save_data_stat||do_stat_always) {
Mdouble fdott = forcet.GetLength();
Mdouble deltat_norm = -deltat.GetLength();
Vec3D centre = 0.5 * (PI->get_Position() + PJ->get_Position());
if (!PI->isFixed())
{
if(save_data_stat||do_stat_always) gather_statistics_collision(PJreal->get_Index(),PI->get_Index(), centre, deltan, deltat_norm,fdotn,fdott,-normal,-(fdott?forcet/fdott:forcet));
if(save_data_fstat) fstat_file
<< t << " "
<< PJreal->get_Index() << " "
<< PI->get_Index() << " "
<< centre << " "
<< deltan << " "
<< deltat_norm << " "
<< fdotn << " "
<< fdott << " "
<< -normal << " "
<< -(fdott?forcet/fdott:forcet) << endl;
}
if (!PJreal->isFixed())
{
if(save_data_stat||do_stat_always) gather_statistics_collision(PI->get_Index(),PJreal->get_Index(), centre, deltan, deltat_norm,fdotn,fdott,normal,(fdott?forcet/fdott:forcet));
if(save_data_fstat){ fstat_file
<< t << " "
<< PI->get_Index() << " "
<< PJreal->get_Index() << " "
<< centre << " "
<< deltan << " "
<< deltat_norm << " "
<< fdotn << " "
<< fdott << " "
<< normal << " "
<< (fdott?forcet/fdott:forcet) << endl;
}
}
}
} // end if particle i and j are overlapping
}
| void MD::compute_walls | ( | Particle * | PI | ) | [protected, virtual] |
This is were the walls are.
This is were the walls are implemented - normals are outward normals.
References Particle::add_Force(), Particle::add_Torque(), CSpecies::disp, CSpecies::disprolling, CSpecies::dispt, CSpecies::disptorsion, do_stat_always, dt, ene_ela, STD_save::fstat_file, gather_statistics_collision(), Particle::get_AngularVelocity(), get_Hertzian(), Particle::get_Index(), Particle::get_IndSpecies(), get_MixedSpecies(), get_NWall(), Particle::get_PeriodicFromParticle(), Particle::get_Position(), Particle::get_Radius(), Particle::get_TangentialSprings(), Particle::get_Torque(), Particle::get_Velocity(), Vec3D::GetLength(), Vec3D::GetLength2, Particle::isFixed(), CSpecies::k, CSpecies::krolling, CSpecies::kt, CSpecies::ktorsion, CSpecies::mu, CSpecies::murolling, CSpecies::mus, CSpecies::musrolling, CSpecies::mustorsion, CSpecies::mutorsion, save_data_ene, save_data_fstat, save_data_stat, CTangentialSprings::select_wall_spring(), Particle::set_Torque(), Vec3D::set_zero(), Species, sqr, t, and Walls.
Referenced by compute_external_forces().
{
//No need to compute interactions between periodic particle images and walls
if(pI->get_PeriodicFromParticle()!=NULL)
return;
for (int i=0; i<get_NWall(); i++)
{
// note: normal points away from particle i
Vec3D normal;
Mdouble dist;
bool touch = Walls[i].get_distance_and_normal(*pI, dist, normal);
//If the wall is being touched (I think :Ant)
if (touch)
{
CSpecies* pSpecies = (pI->get_IndSpecies()==Walls[i].indSpecies)?&Species[pI->get_IndSpecies()]:get_MixedSpecies(pI->get_IndSpecies(),Walls[i].indSpecies);
Mdouble deltan = pI->get_Radius()-dist;
// Compute the relative velocity vector v_ij
Vec3D vrel;
if (!pSpecies->mu) {
vrel = pI->get_Velocity() - Walls[i].get_velocity();
} else {
vrel = pI->get_Velocity() - Walls[i].get_velocity() - Cross(normal, pI->get_AngularVelocity()) * dist;
}
// Compute the projection of vrel onto the normal (can be negative)
Mdouble vdotn = -Dot(vrel, normal);
// Compute normal force on particle i due to contact
Mdouble fdotn = pSpecies->k * deltan - pSpecies->disp * vdotn;
Vec3D force = normal * (-fdotn);
//If tangential forces are present
CTangentialSpring* TangentialSpring = NULL;
Vec3D forcet; forcet.set_zero();
Vec3D forcerolling; forcerolling.set_zero();
Vec3D forcetorsion; forcetorsion.set_zero();
if (pSpecies->mu || pSpecies->murolling || pSpecies->mutorsion) {
//call tangential spring
if (pSpecies->kt || pSpecies->krolling || pSpecies->ktorsion)
TangentialSpring = pI->get_TangentialSprings().select_wall_spring(i, t, dt);
//Compute norm of normal force
Mdouble norm_fn = abs(fdotn);
//calculate sliding friction
if (pSpecies->mu) {
//Compute the tangential component of vrel
Vec3D vrelt=vrel+normal*vdotn;
//Compute norm of vrelt
Mdouble vdott=vrelt.GetLength();
if (pSpecies->kt) {
Vec3D* delta = &(TangentialSpring->delta);
//Integrate the spring
(*delta) += vrelt * dt; //no correction because the normal direction is constant
//Calculate test force including viscous force
forcet = (-pSpecies->dispt) * vrelt - pSpecies->kt * (*delta);
Mdouble forcet2 = forcet.GetLength2();
//tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking),
//but the force is limited by Coulomb friction (sliding):
//f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
if( forcet2 <= sqr(pSpecies->mus*norm_fn) ) {
//sticking++;
TangentialSpring->sliding=false;
} else {
//sliding++;
TangentialSpring->sliding=true;
Mdouble norm_forcet = sqrt(forcet2);
forcet *= pSpecies->mu * norm_fn / norm_forcet;
(*delta) = (forcet + pSpecies->dispt * vrelt)/(-pSpecies->kt);
//~ (*delta) = forcet/(-pSpecies->kt);
}
//Add tangential force to total force
force += forcet;
} else { //if no tangential spring
//tangential forces are modelled by a damper of viscosity dispt (sticking),
//but the force is limited by Coulomb friction (sliding):
//f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
if (vdott*pSpecies->dispt <= pSpecies->mus*norm_fn) { //sticking++;
forcet = -pSpecies->dispt * vrelt;
} else { //sliding++;
//set force to Coulomb limit
forcet = -(pSpecies->mu * norm_fn / vdott) * vrelt;
}
//Add tangential force to total force
force += forcet;
}
}
//calculate rolling friction
if (pSpecies->murolling) {
//From Luding2008, objective rolling velocity (eq 15) w/o 2.0!
Mdouble reducedRadiusI = pI->get_Radius() - .5 * deltan;
Mdouble reducedRadiusIJ = reducedRadiusI;
Vec3D vrolling=-reducedRadiusIJ * Cross(normal, pI->get_AngularVelocity());
if (pSpecies->krolling) {
Vec3D* RollingSpring = &(TangentialSpring->RollingSpring);
//Integrate the spring
(*RollingSpring) += vrolling * dt;// - Dot(*RollingSpring,normal)*normal;
//Calculate test force including viscous force
forcerolling = (-pSpecies->disprolling) * vrolling - pSpecies->krolling * (*RollingSpring);
Mdouble forcerolling2 = forcerolling.GetLength2();
//tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking),
//but the force is limited by Coulomb friction (sliding):
//f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
if( forcerolling2 <= sqr(pSpecies->musrolling*norm_fn) ) {
//sticking++;
} else {
//sliding++;
forcerolling *= pSpecies->murolling * norm_fn / sqrt(forcerolling2);
(*RollingSpring) = (forcerolling + pSpecies->disprolling * vrolling)/(-pSpecies->krolling);
}
//Add tangential force to torque
Vec3D Torque = reducedRadiusIJ * Cross(normal, forcerolling);
pI->add_Torque(Torque);
}
}
//calculate torsive friction
if (pSpecies->mutorsion) {
//From Luding2008, spin velocity (eq 16) w/o 2.0!
Mdouble RadiusIJ = pI->get_Radius();
Vec3D vtorsion=RadiusIJ*Dot(normal,pI->get_AngularVelocity())*normal;
if (pSpecies->ktorsion) {
//~ cout << "Error; not yet implemented" << endl;
//~ exit(-1);
Vec3D* TorsionSpring = &(TangentialSpring->TorsionSpring);
//Integrate the spring
(*TorsionSpring) = Dot((*TorsionSpring) + vtorsion * dt,normal)*normal;
//Calculate test force including viscous force
forcetorsion = (-pSpecies->disptorsion) * vtorsion - pSpecies->ktorsion * (*TorsionSpring);
Mdouble forcetorsion2 = forcetorsion.GetLength2();
//tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking),
//but the force is limited by Coulomb friction (sliding):
//f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
if( forcetorsion2 <= sqr(pSpecies->mustorsion*norm_fn) ) {
//sticking++;
} else {
//sliding++;
cout << "sliding " << endl;
forcetorsion *= pSpecies->mutorsion * norm_fn / sqrt(forcetorsion2);
(*TorsionSpring) = (forcetorsion + pSpecies->disptorsion * vtorsion)/(-pSpecies->ktorsion);
}
//Add tangential force to torque
Vec3D torque = RadiusIJ * forcetorsion;
pI->add_Torque(torque);
}
}
} //end if tangential forces
//Make force Hertzian (note: deltan is normalized by the normal distanct of two particles in contact, as in Silbert)
if (get_Hertzian()) force *= sqrt(deltan/(2.0*pI->get_Radius()));
//Add forces to total force
pI->add_Force(force);
if (pSpecies->mu) pI->set_Torque(pI->get_Torque()+Cross(normal, force) * dist);
// Add torque due to tangential forces: t = Cross(l,f), l=dist*Wall.normal
if (save_data_ene) {
ene_ela += 0.5 * (pSpecies->k * sqr(deltan) + (TangentialSpring?
(pSpecies->kt * TangentialSpring->delta.GetLength2()
+pSpecies->krolling * TangentialSpring->RollingSpring.GetLength2()
+pSpecies->ktorsion * TangentialSpring->TorsionSpring.GetLength2()):0.0));
}
if (save_data_fstat||save_data_stat||do_stat_always)
{
Mdouble fdott = forcet.GetLength();
Mdouble deltat_norm = TangentialSpring?(-TangentialSpring->delta.GetLength()):0.0;
if (!pI->isFixed())
{
if(save_data_stat||do_stat_always) gather_statistics_collision(pI->get_Index(),-(i+1), pI->get_Position() + normal*(pI->get_Radius()-deltan), deltan, deltat_norm,fdotn,fdott,normal,-(fdott?forcet/fdott:forcet));
if(save_data_fstat) fstat_file
<< t << " "
<< pI->get_Index() << " "
<< -(i+1) << " "
<< pI->get_Position() + normal*(pI->get_Radius()-deltan) << " "
<< deltan << " "
<< deltat_norm << " "
<< fdotn << " "
<< fdott << " "
<< normal << " "
<< -(fdott?forcet/fdott:forcet) << endl;
}
}
} // end if particle i touches the wall
}//end if for wall i
}
| void MD::constructor | ( | ) |
A public constructor, which sets defaults so the problem can be solve off the shelf.
todo{to incorporate changes in icc a make clean is required. Why?}
This is where the constructor is defines setup a basic two dimensional problem which can be solved off the shelf///
Reimplemented from STD_save.
Reimplemented in HGRID_base, ChuteWithHopper, Chute, HGRID_3D, HGRID_2D, StatisticsVector< T >, and ChuteBottom.
References data_FixedParticles, dim, dt, format, gravity, Hertzian, particleHandler, PeriodicCreated, STD_save::problem_name, random, save_restart_data_counter, set_dim_particle(), set_do_stat_always(), set_k(), set_NWall(), set_NWallPeriodic(), RNG::set_RandomSeed(), set_restarted(), set_rho(), set_save_count_all(), ParticleHandler::set_StorageCapacity(), Species, tmax, Walls, xballs_additional_arguments, xballs_cmode, xballs_scale, xballs_vscale, xmax, xmin, ymax, ymin, zmax, and zmin.
{
Species.resize(1);
set_restarted(false);
//This sets the maximum number of particles
particleHandler.set_StorageCapacity(2);
dim = 2;
format = 0;
//These are the particle parameters like dissipation etc..
set_k(1e4);
set_rho(2000);
//\todo{Why shouldn't we set dim_particle (which defines the mass calculations) ALWAYS to three, even for 2D problems?}
set_dim_particle(2);
Hertzian=false;
// Set gravitational acceleration
gravity = Vec3D(0.0, -9.8, 0.0);
//This is the parameter of the numerical part
dt=0e-5; // if dt is not user-specified, this is set in actions_before_time_loop()
set_save_count_all(20);
set_do_stat_always(false);
tmax=1.0;
//This sets the default xballs domain
xmin=0.0;
xmax=0.01;
ymin=0.0;
ymax=0.01;
zmin=0.0;
zmax=0.0;
//This set the default position of walls
set_NWall(4);
Walls[0].set(Vec3D(-1,0,0), -xmin);
Walls[1].set(Vec3D( 1,0,0), xmax);
Walls[2].set(Vec3D(0,-1,0), -ymin);
Walls[3].set(Vec3D(0, 1,0), ymax);
// No periodic Walls
set_NWallPeriodic(0);
problem_name << "Problem_1";
data_FixedParticles = 0;
//Default mode is energy with no scale of the vectors
xballs_cmode=0;
xballs_vscale=-1;
xballs_scale=-1;
xballs_additional_arguments = "";
save_restart_data_counter = 0;
PeriodicCreated=0;
//The default random seed is 0
random.set_RandomSeed(0);
#ifdef DEBUG_OUTPUT
cerr << "MD problem constructor finished " << endl;
#endif
}
| virtual bool MD::continue_solve | ( | ) | [inline, protected, virtual] |
Referenced by solve().
{return true;}
| virtual void MD::cout_time | ( | ) | [inline, protected, virtual] |
| virtual void MD::create_xballs_script | ( | ) | [virtual] |
This creates a scipt which can be used to load the xballs problem to display the data just generated.
Referenced by solve().
| void MD::do_integration_after_force_computation | ( | Particle * | pI | ) | [protected, virtual] |
This is were the integration is done.
This is were the integration is done, at the moment it is Velocity Verlet integration.
References Particle::accelerate(), Particle::angularAccelerate(), dt, Particle::get_AngularVelocity(), Particle::get_Force(), get_HGRID_UpdateEachTimeStep(), Particle::get_InvInertia(), Particle::get_InvMass(), get_NWallPeriodic(), Particle::get_Position(), Particle::get_Torque(), HGRID_RemoveParticleFromHgrid(), HGRID_UpdateParticleInHgrid(), Particle::rotate(), rotation, Particle::set_Position(), Particle::set_Velocity(), and WallsPeriodic.
Referenced by solve().
{
#ifdef USE_SIMPLE_VERLET_INTEGRATION
static Vec3D xtemp, atemp;
xtemp=pI->get_Position();
atemp=pI->get_Force()*pI->get_InvMass();
pI->set_Position(xtemp*2.0-pI->PrevPosition+atemp*dt*dt);
pI->set_Velocity((pI->get_Position()-pI->PrevPosition)/(2*dt)+atemp*dt);
pI->PrevPosition=xtemp;
if (rotation) {
pI->angularAccelerate(pI->get_Torque()*pI->get_InvInertia()*dt);
pI->rotate(pI->get_AngularVelocity()*dt);
}
for (int j=0; j<get_NWallPeriodic(); j++)
{
if (WallsPeriodic[j].distance(pI)<0)
{
WallsPeriodic[j].shift_position(pI->get_Position());
if(!get_HGRID_UpdateEachTimeStep())
{
HGRID_RemoveParticleFromHgrid(iP);
HGRID_UpdateParticleInHgrid(iP);
}
}
}
#else //use velocity verlet
pI->accelerate(pI->get_Force()*pI->get_InvMass()*0.5*dt);
if (rotation) pI->get_AngularVelocity() += pI->get_Torque()*pI->get_InvInertia()*0.5*dt;
#endif
}
| void MD::do_integration_before_force_computation | ( | Particle * | pI | ) | [protected, virtual] |
This is were the integration is done.
This is were the integration is done, at the moment it is Velocity Verlet integration.
References Particle::accelerate(), Particle::angularAccelerate(), dt, Particle::get_AngularVelocity(), Particle::get_Force(), get_HGRID_UpdateEachTimeStep(), Particle::get_InvInertia(), Particle::get_InvMass(), get_NWallPeriodic(), Particle::get_Position(), Particle::get_Torque(), Particle::get_Velocity(), Vec3D::GetLength(), HGRID_RemoveParticleFromHgrid(), HGRID_update_move(), HGRID_UpdateParticleInHgrid(), Particle::move(), Particle::rotate(), rotation, and WallsPeriodic.
Referenced by solve().
{
#ifdef USE_SIMPLE_VERLET_INTEGRATION
#else //use velocity verlet
iP->accelerate(iP->get_Force()*iP->get_InvMass()*0.5*dt);
iP->move(iP->get_Velocity()*dt);
HGRID_update_move(iP,iP->get_Velocity().GetLength()*dt);
if (rotation)
{
iP->angularAccelerate(iP->get_Torque()*iP->get_InvInertia()*0.5*dt);
iP->rotate(iP->get_AngularVelocity()*dt);
}
// This shifts particles that moved through periodic walls
for (int j=0; j<get_NWallPeriodic(); j++)
{
if (WallsPeriodic[j].distance(*iP)<0)
{
WallsPeriodic[j].shift_position(iP->get_Position());
if(!get_HGRID_UpdateEachTimeStep())
{
HGRID_RemoveParticleFromHgrid(iP);
HGRID_UpdateParticleInHgrid(iP);
}
}
}
#endif
}
| bool MD::find_next_data_file | ( | Mdouble | tmin, |
| bool | verbose = true |
||
| ) |
References STD_save::data_file, STD_save::get_data_filename(), STD_save::get_file_counter(), STD_save::get_options_data(), STD_save::increase_counter_data(), STD_save::set_file_counter(), and t.
{
if (get_options_data()==2) {
while(true) {
increase_counter_data(fstream::in);
//check if file exists and contains data
int N;
data_file >> N;
if (data_file.eof()||data_file.peek()==-1) {
cout << "file " << get_data_filename() << " not found" << endl;
return false;
}
//check if tmin condition is satisfied
Mdouble t;
data_file >> t;
if (t>tmin) {
set_file_counter(get_file_counter()-1);
return true;
}
if (verbose) cout <<"Jumping file counter: "<<get_file_counter() << endl;
}
} else return true;
//length = is.tellg();
//is.seekg (0, ios::end);
}
| virtual void MD::finish_statistics | ( | ) | [inline, protected, virtual] |
| void MD::fstat_header | ( | ) | [protected, virtual] |
References STD_save::fstat_file, get_LargestParticle(), Particle::get_Radius(), get_SmallestParticle(), get_t(), get_xmax(), get_xmin(), get_ymax(), get_ymin(), get_zmax(), and get_zmin().
Referenced by solve().
{
// line #1: time, volume fraction
// line #2: wall box: wx0, wy0, wz0, wx1, wy1, wz1
// line #3: radii-min-max & moments: rad_min, rad_max, r1, r2, r3, r4
fstat_file << "#"
<< " " << get_t()
<< " " << 0
<< endl;
fstat_file << "#"
<< " " << get_xmin()
<< " " << get_ymin()
<< " " << get_zmin()
<< " " << get_xmax()
<< " " << get_ymax()
<< " " << get_zmax()
<< endl;
fstat_file << "#"
<< " " << get_SmallestParticle()->get_Radius()
<< " " << get_LargestParticle()->get_Radius()
<< " " << 0
<< " " << 0
<< " " << 0
<< " " << 0
<< endl;
}
| virtual void MD::gather_statistics_collision | ( | int index1 | UNUSED, |
| int index2 | UNUSED, | ||
| Vec3D Contact | UNUSED, | ||
| Mdouble delta | UNUSED, | ||
| Mdouble ctheta | UNUSED, | ||
| Mdouble fdotn | UNUSED, | ||
| Mdouble fdott | UNUSED, | ||
| Vec3D P1_P2_normal_ | UNUSED, | ||
| Vec3D P1_P2_tangential | UNUSED | ||
| ) | [inline, protected, virtual] |
Reimplemented in StatisticsVector< T >.
Referenced by compute_internal_forces(), compute_plastic_internal_forces(), and compute_walls().
{};
| bool MD::get_append | ( | ) | [inline] |
| Mdouble MD::get_collision_time | ( | Mdouble | mass, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Calculates collision time for two copies of a particle of given disp, k, mass.
{return Species[indSpecies].get_collision_time(mass);}
| Mdouble MD::get_depth | ( | unsigned int | indSpecies = 0 | ) | [inline] |
{return Species[indSpecies].get_depth();}
| int MD::get_dim | ( | ) | [inline] |
Allows the dimension of the simulation to be accessed.
Referenced by StatisticsVector< T >::set_nz().
{return dim;}
| int MD::get_dim_particle | ( | unsigned int | indSpecies = 0 | ) | [inline] |
Allows the dimension of the particle (f.e. for mass) to be accessed.
{return Species[indSpecies].get_dim_particle();}
| Mdouble MD::get_disp | ( | unsigned int | indSpecies = 0 | ) | [inline] |
Allows the normal dissipation to be accessed.
Referenced by readNextArgument().
{return Species[indSpecies].get_dissipation();}
| Mdouble MD::get_disprolling | ( | unsigned int | indSpecies = 0 | ) | [inline] |
Allows the tangential viscosity to be accessed.
{return Species[indSpecies].get_disprolling();}
| Mdouble MD::get_dispt | ( | unsigned int | indSpecies = 0 | ) | [inline] |
Allows the tangential viscosity to be accessed.
Referenced by readNextArgument().
{return Species[indSpecies].get_dispt();}
| Mdouble MD::get_disptorsion | ( | unsigned int | indSpecies = 0 | ) | [inline] |
Allows the tangential viscosity to be accessed.
{return Species[indSpecies].get_disptorsion();}
| Mdouble MD::get_dissipation | ( | unsigned int | indSpecies = 0 | ) | [inline] |
Allows the normal dissipation to be accessed.
Referenced by Chute::set_collision_time_and_restitution_coefficient().
{return Species[indSpecies].get_dissipation();}
| int MD::get_do_stat_always | ( | ) | [inline] |
{return do_stat_always;}
| Mdouble MD::get_dt | ( | ) | [inline] |
Allows the time step dt to be accessed.
Referenced by StatisticsVector< T >::check_current_time_for_statistics(), load_par_ini_file(), ChuteBottom::make_rough_bottom(), and readNextArgument().
{return dt;}
| Mdouble MD::get_ene_ela | ( | ) | [inline] |
Gets ene_ela.
{return ene_ela;}
| Vec3D MD::get_gravity | ( | ) | [inline] |
Allows the gravitational acceleration to be accessed.
Referenced by Chute::set_ChuteAngle().
{return gravity;}
| bool MD::get_Hertzian | ( | ) | [inline] |
Referenced by compute_internal_forces(), and compute_walls().
{return Hertzian;}
| virtual bool MD::get_HGRID_UpdateEachTimeStep | ( | ) | [inline, protected, virtual] |
Reimplemented in HGRID_base.
Referenced by do_integration_after_force_computation(), and do_integration_before_force_computation().
{return true;};
Allows the spring constant to be accessed.
Referenced by readNextArgument().
{return Species[indSpecies].get_k();}
| Mdouble MD::get_k1 | ( | unsigned int | indSpecies = 0 | ) | [inline] |
Allows the plastic constants to be accessed.
{return Species[indSpecies].get_k1();}
| Mdouble MD::get_k2max | ( | unsigned int | indSpecies = 0 | ) | [inline] |
{return Species[indSpecies].get_k2max();}
| Mdouble MD::get_kc | ( | unsigned int | indSpecies = 0 | ) | [inline] |
{return Species[indSpecies].get_kc();}
| Mdouble MD::get_krolling | ( | int | indSpecies = 0 | ) | [inline] |
Allows the spring constant to be accessed.
{return Species[indSpecies].get_krolling();}
| Mdouble MD::get_kt | ( | int | indSpecies = 0 | ) | [inline] |
Allows the spring constant to be accessed.
Referenced by readNextArgument().
{return Species[indSpecies].get_kt();}
| Mdouble MD::get_ktorsion | ( | int | indSpecies = 0 | ) | [inline] |
Allows the spring constant to be accessed.
{return Species[indSpecies].get_ktorsion();}
| virtual Particle* MD::get_LargestParticle | ( | ) | [inline, virtual] |
Reimplemented in Chute.
Referenced by HGRID_base::fix_hgrid(), fstat_header(), and solve().
{return particleHandler.get_LargestParticle();}
| Mdouble MD::get_Mass_from_Radius | ( | Mdouble | radius, |
| int | indSpecies = 0 |
||
| ) | [inline] |
References Particle::compute_particle_mass(), Particle::get_Mass(), Particle::set_IndSpecies(), and Particle::set_Radius().
{
Particle P;
P.set_Radius(radius);
P.set_IndSpecies(indSpecies);
P.compute_particle_mass(Species);
return P.get_Mass();
}
| Mdouble MD::get_max_radius | ( | ) | [inline] |
| Mdouble MD::get_maximum_velocity | ( | Particle & | P | ) | [inline] |
Calculates the maximum velocity allowed for a collision of two copies of P (for higher velocities particles could pass through each other)
References Particle::get_IndSpecies(), Particle::get_Mass(), and Particle::get_Radius().
{return Species[P.get_IndSpecies()].get_maximum_velocity(P.get_Radius(),P.get_Mass());}
| Mdouble MD::get_maximum_velocity | ( | ) | [inline] |
References Particle::get_IndSpecies(), Particle::get_Mass(), and Particle::get_Radius().
{
Particle *P = get_SmallestParticle();
return P->get_Radius() * sqrt(Species[P->get_IndSpecies()].k/(.5*P->get_Mass()));
}
| CSpecies* MD::get_MixedSpecies | ( | int | i, |
| int | j | ||
| ) | [inline] |
Allows the mixed species to be accessed.
References CSpecies::MixedSpecies.
Referenced by compute_internal_forces(), compute_plastic_internal_forces(), and compute_walls().
| Mdouble MD::get_mu | ( | unsigned int | indSpecies = 0 | ) | [inline] |
Allows the Coulomb friction coefficient to be accessed.
Referenced by readNextArgument(), and solve().
{return Species[indSpecies].get_mu();}
| Mdouble MD::get_murolling | ( | unsigned int | indSpecies = 0 | ) | [inline] |
Allows the Coulomb friction coefficient to be accessed.
{return Species[indSpecies].get_murolling();}
| Mdouble MD::get_mutorsion | ( | unsigned int | indSpecies = 0 | ) | [inline] |
Allows the Coulomb friction coefficient to be accessed.
{return Species[indSpecies].get_mutorsion();}
| int MD::get_NSpecies | ( | ) | [inline] |
Allows the number of Species to be accessed.
{return Species.size();}
| int MD::get_NWall | ( | ) | [inline] |
Allows the number of walls to be accessed.
Referenced by ChuteWithHopper::add_hopper(), compute_walls(), and print().
{return Walls.size();}
| int MD::get_NWallPeriodic | ( | ) | [inline] |
Allows the number of periodic walls to be accessed.
Referenced by Check_and_Duplicate_Periodic_Particle(), do_integration_after_force_computation(), do_integration_before_force_computation(), and print().
{return WallsPeriodic.size();}
| ParticleHandler& MD::get_ParticleHandler | ( | ) | [inline] |
Referenced by ChuteWithHopper::add_hopper(), Chute::add_particle(), Chute::clean_chute(), Chute::cout_time(), Chute::create_bottom(), HGRID_base::fix_hgrid(), Chute::get_LargestParticle(), Chute::get_radius_of_largest_particle(), Chute::get_radius_of_smallest_particle(), Chute::get_SmallestParticle(), HGRID_base::HGRID_actions_before_time_step(), HGRID_base::InitBroadPhase(), Chute::IsInsertable(), ChuteBottom::make_rough_bottom(), HGRID_base::set_HGRID_num_buckets_to_power(), ChuteBottom::setup_particles_initial_conditions(), and solve().
{ return particleHandler;}
| Mdouble MD::get_plastic_dt | ( | Mdouble | mass, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
{return Species[indSpecies].get_plastic_dt(mass);}
| int MD::get_restart_version | ( | ) | [inline] |
Gets restart_version.
Referenced by Chute::read(), and HGRID_base::read().
{return restart_version;}
| bool MD::get_restarted | ( | ) | [inline] |
Gets restarted.
{return restarted;}
| Mdouble MD::get_restitution_coefficient | ( | Mdouble | mass, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Calculates restitution coefficient for two copies of given disp, k, mass.
{return Species[indSpecies].get_restitution_coefficient(mass);}
| Mdouble MD::get_rho | ( | int | indSpecies = 0 | ) | [inline] |
Allows the density to be accessed.
Referenced by Chute::set_collision_time_and_restitution_coefficient().
{return Species[indSpecies].get_rho();}
| int MD::get_save_count | ( | ) | [inline] |
{return get_save_count_data ();}
| int MD::get_save_count_data | ( | ) | [inline] |
{return save_count_data;}
| int MD::get_save_count_ene | ( | ) | [inline] |
{return save_count_ene;}
| int MD::get_save_count_fstat | ( | ) | [inline] |
{return save_count_fstat;}
| int MD::get_save_count_stat | ( | ) | [inline] |
{return save_count_stat;}
| int MD::get_save_data_data | ( | ) | [inline] |
Returns the data counter.
{return save_data_data;}
| int MD::get_save_data_ene | ( | ) | [inline] |
{return save_data_ene;}
| int MD::get_save_data_fstat | ( | ) | [inline] |
{return save_data_fstat;}
| int MD::get_save_data_stat | ( | ) | [inline] |
{return save_data_stat;}
| int MD::get_savecount | ( | ) | [inline] |
Allows the number of time steps between saves to be accessed.
{return get_save_count_data ();}
| virtual Particle* MD::get_SmallestParticle | ( | ) | [inline, virtual] |
Reimplemented in Chute.
Referenced by HGRID_base::fix_hgrid(), fstat_header(), and HGRID_base::InitBroadPhase().
{return particleHandler.get_SmallestParticle();}
| vector<CSpecies>& MD::get_Species | ( | void | ) | [inline] |
Allows the species to be copied.
{return Species;}
| CSpecies* MD::get_Species | ( | int | i | ) | [inline] |
Allows the species to be accessed.
{return &Species[i];}
Access function for the time.
Referenced by StatisticsVector< T >::check_current_time_for_statistics(), Chute::cout_time(), fstat_header(), and output_ene().
{return t;}
| Mdouble MD::get_tmax | ( | ) | [inline] |
Allows the upper time limit to be accessed.
Referenced by Chute::cout_time(), StatisticsVector< T >::get_tmaxStat(), and readNextArgument().
{return tmax;}
| vector<CWall>& MD::get_Walls | ( | void | ) | [inline] |
Allows the walls to be copied.
{return Walls;}
| CWall& MD::get_Walls | ( | int | i | ) | [inline] |
Allows the walls to be accessed.
{return Walls[i];}
| vector<CWallPeriodic>& MD::get_WallsPeriodic | ( | void | ) | [inline] |
Allows the periodic walls to be accessed.
{return WallsPeriodic;}
| CWallPeriodic& MD::get_WallsPeriodic | ( | int | i | ) | [inline] |
Allows the periodic walls to be accessed.
{return WallsPeriodic[i];}
| string MD::get_xballs_additional_arguments | ( | ) | [inline] |
{return xballs_additional_arguments;}
| int MD::get_xballs_cmode | ( | ) | [inline] |
{return xballs_cmode;}
| double MD::get_xballs_scale | ( | ) | [inline] |
{return xballs_scale;}
| double MD::get_xballs_vscale | ( | ) | [inline] |
{return xballs_vscale;}
| Mdouble MD::get_xmax | ( | ) | [inline] |
Get xmax.
Referenced by Chute::clean_chute(), Chute::create_bottom(), fstat_header(), Chute::get_ChuteLength(), ChuteWithHopper::get_ChuteLength(), ChuteWithHopper::get_MaximumVelocityInducedByGravity(), StatisticsVector< T >::get_xmaxStat(), load_par_ini_file(), ChuteWithHopper::set_shift(), and ChuteBottom::setup_particles_initial_conditions().
{return xmax;}
| Mdouble MD::get_xmin | ( | ) | [inline] |
Get xmin.
Referenced by Chute::clean_chute(), Chute::create_bottom(), Chute::create_inflow_particle(), fstat_header(), StatisticsVector< T >::get_xminStat(), load_par_ini_file(), and ChuteBottom::setup_particles_initial_conditions().
{return xmin;}
| Mdouble MD::get_ymax | ( | ) | [inline] |
Gets ymax.
Referenced by ChuteWithHopper::add_hopper(), Chute::create_bottom(), Chute::create_inflow_particle(), ChuteWithHopper::create_inflow_particle(), fstat_header(), Chute::get_ChuteWidth(), StatisticsVector< T >::get_ymaxStat(), load_par_ini_file(), Chute::setup_particles_initial_conditions(), and ChuteBottom::setup_particles_initial_conditions().
{return ymax;}
| Mdouble MD::get_ymin | ( | ) | [inline] |
Gets ymin.
Referenced by ChuteWithHopper::add_hopper(), Chute::create_bottom(), Chute::create_inflow_particle(), ChuteWithHopper::create_inflow_particle(), fstat_header(), StatisticsVector< T >::get_yminStat(), load_par_ini_file(), Chute::setup_particles_initial_conditions(), and ChuteBottom::setup_particles_initial_conditions().
{return ymin;}
| Mdouble MD::get_zmax | ( | ) | [inline] |
Gets zmax.
Referenced by Chute::create_inflow_particle(), fstat_header(), StatisticsVector< T >::get_zmaxStat(), load_par_ini_file(), and ChuteBottom::setup_particles_initial_conditions().
{return zmax;}
| Mdouble MD::get_zmin | ( | ) | [inline] |
Gets zmin.
Referenced by Chute::create_bottom(), Chute::create_inflow_particle(), fstat_header(), StatisticsVector< T >::get_zminStat(), load_par_ini_file(), and ChuteBottom::setup_particles_initial_conditions().
{return zmin;}
| virtual double MD::getInfo | ( | Particle &P | UNUSED | ) | [inline, virtual] |
Allows the user to set what is written into the info column in the data file.
{return 0;}
| virtual void MD::HGRID_actions_after_integration | ( | ) | [inline, protected, virtual] |
| virtual void MD::HGRID_actions_before_integration | ( | ) | [inline, protected, virtual] |
| virtual void MD::HGRID_actions_before_time_loop | ( | ) | [inline, protected, virtual] |
This is actions before the start of the main time loop.
Reimplemented in HGRID_base.
Referenced by solve(), and statistics_from_restart_data().
{};
| virtual void MD::HGRID_actions_before_time_step | ( | ) | [inline, protected, virtual] |
This is action before the time step is started.
Reimplemented in HGRID_base.
Referenced by solve(), and statistics_from_restart_data().
{};
| virtual void MD::HGRID_InsertParticleToHgrid | ( | Particle *obj | UNUSED | ) | [inline, protected, virtual] |
This is action before the time step is started.
Referenced by Check_and_Duplicate_Periodic_Particle().
{};
| virtual void MD::HGRID_RemoveParticleFromHgrid | ( | Particle *obj | UNUSED | ) | [inline, protected, virtual] |
Referenced by do_integration_after_force_computation(), and do_integration_before_force_computation().
{};
| virtual void MD::HGRID_update_move | ( | Particle * | , |
| Mdouble | |||
| ) | [inline, protected, virtual] |
| virtual void MD::HGRID_UpdateParticleInHgrid | ( | Particle *obj | UNUSED | ) | [inline, protected, virtual] |
| void MD::info | ( | ) | [inline, virtual] |
| virtual void MD::InitBroadPhase | ( | ) | [inline, protected, virtual] |
Initialisation of Broad Phase Information (Default no Broad Phase so empty)
Reimplemented in HGRID_base.
{};
| virtual void MD::initialize_statistics | ( | ) | [inline, protected, virtual] |
Functions for statistics.
Reimplemented in StatisticsVector< T >.
Referenced by solve(), and statistics_from_restart_data().
{};
| void MD::initialize_tangential_springs | ( | ) | [protected] |
| bool MD::load_from_data_file | ( | const char * | filename, |
| unsigned int | format = 0 |
||
| ) |
This allows particle data to be reloaded from data files.
This function loads data from the .data file i.e.
you get position, mass and velocity information only See also MD::load_restart_data Can read in format 14 - 8 or format 7 data This code saves in format 8 for 2D and format 14 for 3D so if no extra parameters are specified it will assume this note: many parameters, like density cannot be set using the data file
References STD_save::data_file, STD_save::get_options_data(), STD_save::options_data, read_next_from_data_file(), and STD_save::set_options_data().
Referenced by statistics_from_restart_data().
{
//This opens the file the data will be recalled from
data_file.open( filename, fstream::in);
if (!data_file.is_open()||data_file.bad()) {
cout << "Loading data file " << filename << " failed" << endl;
return false;
}
//options_data is ignored
int options_data = get_options_data();
set_options_data(1);
read_next_from_data_file(format);
set_options_data(options_data);
data_file.close();
cout << "Loaded data file " << filename << endl;
return true;
}
| bool MD::load_par_ini_file | ( | const char * | filename | ) |
allows the user to read par.ini files (useful to read MDCLR files)
References get_dt(), get_xmax(), get_xmin(), get_ymax(), get_ymin(), get_zmax(), get_zmin(), set_dim_particle(), set_disp(), set_dt(), set_gravity(), set_k(), set_NWall(), set_NWallPeriodic(), set_rho(), set_save_count_data(), set_save_count_fstat(), set_savecount(), set_t(), set_tmax(), Walls, and WallsPeriodic.
{
//Opens the par.ini file
fstream file;
file.open( filename, fstream::in);
if (!file.is_open()||file.bad()) {
cout << "Loading par.ini file " << filename << " failed" << endl;
return false;
}
Mdouble doubleValue;
int integerValue;
// inputfile par.ini
// line 1 =============================================================
// Example: 1 1 0
// 1: integer (0|1) switches from non-periodic to periodic
// integer (5|6) does 2D integration only (y-coordinates fixed)
// and switches from non-periodic to periodic
// integer (11) uses a quarter system with circular b.c.
file >> integerValue;
//~ cout << "11" << integerValue << endl;
if (integerValue==0) {
set_NWall(6);
set_NWallPeriodic(0);
Walls[0].set(Vec3D(-1,0,0), -get_xmin());
Walls[1].set(Vec3D( 1,0,0), get_xmax());
Walls[2].set(Vec3D(0,-1,0), -get_ymin());
Walls[3].set(Vec3D(0, 1,0), get_ymax());
Walls[4].set(Vec3D(0,0,-1), -get_zmin());
Walls[5].set(Vec3D(0,0,1), get_zmax());
} else if (integerValue==1) {
set_NWall(0);
set_NWallPeriodic(3);
WallsPeriodic[0].set(Vec3D(1,0,0), get_xmin(), get_xmax());
WallsPeriodic[1].set(Vec3D(0,1,0), get_ymin(), get_ymax());
WallsPeriodic[2].set(Vec3D(0,0,1), get_zmin(), get_zmax());
} else if (integerValue==5) {
set_NWall(4);
set_NWallPeriodic(0);
Walls[0].set(Vec3D(-1,0,0), -get_xmin());
Walls[1].set(Vec3D( 1,0,0), get_xmax());
Walls[2].set(Vec3D(0,-1,0), -get_ymin());
Walls[3].set(Vec3D(0, 1,0), get_ymax());
} else if (integerValue==6) {
set_NWall(0);
set_NWallPeriodic(2);
WallsPeriodic[0].set(Vec3D(1,0,0), get_xmin(), get_xmax());
WallsPeriodic[1].set(Vec3D(0,1,0), get_ymin(), get_ymax());
} else {
cerr << "Error in par.ini: line 1, value 1 is " << integerValue << endl;
exit(-1);
}
// 2: integer (0|1) dont use | use the search pattern for linked cells
file >> integerValue; //ignore
// 3: real - gravity in z direction: positive points upwards
file >> doubleValue;
set_gravity(Vec3D(0.0,0.0,doubleValue));
// line 2 =============================================================
// Example: -10000 .5e-2
// 1: time end of simulation - (negative resets start time to zero
// and uses -time as end-time)
file >> doubleValue;
if (doubleValue<0) set_t(0);
set_tmax(abs(doubleValue));
// 2: time-step of simulation
file >> doubleValue;
set_dt(doubleValue);
// line 3 =============================================================
// Example: 1e-1 100
file >> doubleValue;
if (doubleValue>=0) {
// 1: time-step for output on time-series protocoll file -> "ene"
int savecount = round(doubleValue/get_dt());
set_savecount(savecount);
// 2: time-step for output on film (coordinate) file -> "c3d"
// (fstat-output is coupled to c3d-output time-step)
file >> doubleValue;
savecount = round(doubleValue/get_dt());
set_save_count_data(savecount);
set_save_count_fstat(savecount);
} else {
// or: ---------------------------------------------------------------
// 1: negative number is multiplied to the previous log-output-time
// 2: requires initial log-output time
// 3: negative number is multiplied to the previous film-output-time
// 4: requires initial film-output time
cerr << "Error in par.ini: line 3, value 1 is " << doubleValue << endl;
exit(-1);
}
// line 4 =============================================================
// Example: 2000 1e5 1e3 1e2
// 1: particle density (mass=4/3*constants::pi*density*rad^3)
file >> doubleValue;
set_dim_particle(3);
set_rho(doubleValue);
// 2: linear spring constant
file >> doubleValue;
set_k(doubleValue);
// 3: linear dashpot constant
file >> doubleValue;
set_disp(doubleValue);
// 4: background damping dashpot constant
file >> doubleValue;
if (doubleValue) cerr << "Warning in par.ini: ignored background damping " << doubleValue << endl;
// line 5 =============================================================
// Example: 0 0
// 1: growth rate: d(radius) = xgrow * dt
file >> doubleValue;
if (doubleValue) cerr << "Warning in par.ini: ignored growth rate " << doubleValue << endl;
// 2: target volume_fraction
file >> doubleValue;
if (doubleValue) cerr << "Warning in par.ini: ignored target volume_fraction " << doubleValue << endl;
file.close();
cout << "Loaded par.ini file " << filename << endl;
return true;
}
| int MD::load_restart_data | ( | ) |
Loads all MD data.
This function loads all MD data - See also MD::save_restart_data, MD::load_from_data_file This function return 1 if sucessful else it returns -1.
References STD_save::problem_name.
Referenced by readNextArgument(), and statistics_from_restart_data().
{
stringstream restart_filename;
restart_filename.str("");
restart_filename << problem_name.str().c_str() << ".restart";
return load_restart_data(restart_filename.str());
}
| int MD::load_restart_data | ( | string | filename | ) |
References read(), and set_restarted().
{
ifstream restart_data(filename.c_str());
if( restart_data.good() )
{
read(restart_data);
restart_data.close();
set_restarted(true);
return(1);
} else {
cerr << "restart_data " << filename << " could not be loaded from " << filename << endl;
exit(-1);
}
}
| void MD::output_ene | ( | ) | [protected, virtual] |
This function outputs statistical data - The default is to compute the rotational kinetic engergy, linear kinetic energy, and the centre of mass.
todo{Why is there a +6 here?}
References ParticleHandler::begin(), ParticleHandler::end(), ene_ela, STD_save::ene_file, get_t(), gravity, and particleHandler.
Referenced by solve(), and statistics_from_restart_data().
{
Mdouble ene_kin = 0, ene_rot = 0, ene_gra = 0, mass_sum= 0, x_masslength=0, y_masslength=0, z_masslength=0;
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) if (!(*it)->isFixed())
{
ene_kin += .5 * (*it)->get_Mass() * (*it)->get_Velocity().GetLength2();
ene_rot += .5 * (*it)->get_Inertia() * (*it)->get_AngularVelocity().GetLength2();
ene_gra -= (*it)->get_Mass() * Dot(gravity,(*it)->get_Position());
mass_sum +=(*it)->get_Mass();
x_masslength +=(*it)->get_Mass()*(*it)->get_Position().X;
y_masslength +=(*it)->get_Mass()*(*it)->get_Position().Y;
z_masslength +=(*it)->get_Mass()*(*it)->get_Position().Z;
} //end for loop over Particles
static int width = ene_file.precision() + 6;
ene_file << setw(width) << get_t()
<< " " << setw(width) << ene_gra
<< " " << setw(width) << ene_kin
<< " " << setw(width) << ene_rot
<< " " << setw(width) << ene_ela
<< " " << setw(width) << (mass_sum?x_masslength/mass_sum:NAN)
<< " " << setw(width) << (mass_sum?y_masslength/mass_sum:NAN)
<< " " << setw(width) << (mass_sum?z_masslength/mass_sum:NAN)
<< endl;
ene_ela = 0;
//sliding = sticking = 0;
}
| virtual void MD::output_statistics | ( | ) | [inline, protected, virtual] |
| void MD::output_xballs_data | ( | ) | [protected, virtual] |
Output xball data for Particle i.
This function outputs the location and velocity of the particle in a format the xballs progream can read.
References STD_save::data_file, dim, format, ParticleHandler::get_NumberOfParticles(), output_xballs_data_particle(), particleHandler, t, xmax, xmin, ymax, ymin, zmax, and zmin.
Referenced by solve().
{
//Set the correct formation based of dimension if the formation is not specified by the user
if (format ==0) {
switch (dim) {
case 1:
case 2:
format=8;
break;
case 3:
format=14;
break;
}
}
// This outputs the location of walls and how many particles there are to file this is required by the xballs plotting
if (format!=14) // dim = 1 or 2
data_file << particleHandler.get_NumberOfParticles() << " " <<t << " "
<< xmin << " " << ymin << " "
<< xmax << " " << ymax << " " << endl;
else //dim==3
data_file << particleHandler.get_NumberOfParticles() << " " <<t << " "
<< xmin << " " << ymin << " " << zmin << " "
<< xmax << " " << ymax << " " << zmax << " " << endl;
// This outputs the particle data
for (unsigned int i = 0;i<particleHandler.get_NumberOfParticles();i++)
output_xballs_data_particle(i);
#ifdef DEBUG_OUTPUT
cerr << "Have output the properties of the problem to disk " << endl;
#endif
}
| virtual void MD::output_xballs_data_particle | ( | int | i | ) | [protected, virtual] |
Referenced by output_xballs_data().
| void MD::print | ( | std::ostream & | os, |
| bool | print_all = false |
||
| ) | [virtual] |
Outputs MD.
Reimplemented in HGRID_base, and Chute.
References dim, dt, ParticleHandler::get_NumberOfParticles(), get_NWall(), get_NWallPeriodic(), ParticleHandler::get_Particle(), ParticleHandler::get_StorageCapacity(), gravity, STD_save::options_data, STD_save::options_ene, STD_save::options_fstat, STD_save::options_restart, particleHandler, Particle::print(), STD_save::problem_name, save_count_data, save_count_ene, save_count_fstat, save_count_stat, Species, t, tmax, Walls, WallsPeriodic, xmax, xmin, ymax, ymin, zmax, and zmin.
Referenced by StatisticsVector< T >::statistics_from_fstat_and_data().
{
os << "problem_name:" << problem_name.str().c_str() << endl;
os << " t:" << t << " dt:" << dt << ", tmax:" << tmax << ", save_count_data:" << save_count_data << ", save_count_ene:" << save_count_ene << ", save_count_stat:" << save_count_stat << ", save_count_fstat:" << save_count_fstat << endl;
os << " x:[" << xmin << "," << xmax << "], y:[" << ymin << "," << ymax << "], z:[" << zmin << "," << zmax << "]" << endl;
os << " dim:" << dim << ", gravity:" << gravity << endl;
if (Species.size()==1) {
os << " "; Species[0].print(os); os << endl;
} else {
os << " Species: size:" << Species.size() << endl;
for (unsigned int i=0; i<Species.size(); i++) {
os << " "; Species[i].print(os); os << endl;
for (unsigned int j=0; j<Species[i].MixedSpecies.size(); j++) { os << " "; Species[i].MixedSpecies[j].print(os); os << endl; }
}
}
os << " options_fstat=" << options_fstat
<< " , options_data=" << options_data
<< " , options_ene=" << options_ene
<< " , options_restart=" << options_restart << endl;
os << " Particles: N:" << particleHandler.get_NumberOfParticles() << ", Nmax:" << particleHandler.get_StorageCapacity() << endl;
if (print_all) for (unsigned int i=0; i<particleHandler.get_NumberOfParticles(); i++) { os << " "; particleHandler.get_Particle(i)->print(os); os << endl; }
os << " Walls: NWall:" << get_NWall() << ", NWallPeriodic:" << get_NWallPeriodic() << endl;
for (int i=0; i<get_NWall(); i++) { os << " "; Walls[i].print(os); os << endl; }
for (int i=0; i<get_NWallPeriodic(); i++) { os << " "; WallsPeriodic[i].print(os,print_all); os << endl; }
}
| virtual void MD::process_statistics | ( | bool usethese | UNUSED | ) | [inline, protected, virtual] |
| void MD::read | ( | std::istream & | is | ) | [virtual] |
Reads all MD data.
Reimplemented in ChuteWithHopper, HGRID_base, and Chute.
References Particle::compute_particle_mass(), ParticleHandler::copyAndAddParticle(), data_FixedParticles, dim, dt, gravity, STD_save::options_data, STD_save::options_ene, STD_save::options_fstat, STD_save::options_restart, particleHandler, STD_save::problem_name, read_v1(), read_v2(), restart_version, save_count_data, save_count_ene, save_count_fstat, save_count_stat, Species, t, tmax, Walls, WallsPeriodic, xmax, xmin, ymax, ymin, zmax, and zmin.
Referenced by load_restart_data().
{
string dummy;
is >> dummy;
if (strcmp(dummy.c_str(),"restart_version"))
{
dim = atoi(dummy.c_str());
restart_version = 1;
read_v1(is);
}
else
{
is >> restart_version
>> dummy >> dummy;
problem_name.str(dummy);
if(restart_version==2)
{
read_v2(is);
}
else
{
cout << "version 3" << endl;
is >> dummy >> xmin
>> dummy >> xmax
>> dummy >> ymin
>> dummy >> ymax
>> dummy >> zmin
>> dummy >> zmax
>> dummy >> dt
>> dummy >> t
>> dummy >> tmax
>> dummy >> save_count_data
>> dummy >> save_count_ene
>> dummy >> save_count_stat
>> dummy >> save_count_fstat
>> dummy >> dim
>> dummy >> gravity
>> dummy >> options_fstat
>> dummy >> options_data
>> dummy >> options_ene
>> dummy;
//this is optional to allow restart files with and without options_restart
if (!strcmp(dummy.c_str(),"options_restart")) is >> options_restart >> dummy;
unsigned int N;
is >> N;
Species.resize(N);
//please don't change this code to iterators; the mixed species requires indexing
//old code: for (vector<CSpecies>::iterator it = Species.begin(); it!=Species.end(); ++it) it->read(is);
for (unsigned int i=0; i<N; i++) {
Species[i].read(is);
Species[i].MixedSpecies.resize(i);
for (unsigned int j=0; j<i; j++) {
Species[i].MixedSpecies[j].read(is);
}
}
Species.resize(N);
is >> dummy >> N;
Walls.resize(N);
for (vector<CWall>::iterator it = Walls.begin(); it!=Walls.end(); ++it) it->read(is);
is >> dummy >> N;
WallsPeriodic.resize(N);
for (vector<CWallPeriodic>::iterator it = WallsPeriodic.begin(); it!=WallsPeriodic.end(); ++it) it->read(is);
is >> dummy >> N;
Particle p0;
//for each line of particle data, the line is first completely read. This is to read restart datafiles both with and without tangential spring data.
string line_string;
getline(is,line_string);
for (unsigned int i=0;i<N;i++) {
getline(is,line_string);
stringstream line (stringstream::in | stringstream::out);
line << line_string;
line >> p0;
p0.compute_particle_mass(Species);
particleHandler.copyAndAddParticle(p0);
}
//After the lines storing particle data, an optional variable 'FixedParticles' can be read in, which fixes the first 'FixedParticles' particles
if (is.peek()==70) is >> dummy >> data_FixedParticles;
}
}
}
| int MD::read_dim_from_data_file | ( | ) |
References STD_save::data_file, STD_save::data_filename, t, xmax, xmin, ymax, ymin, zmax, and zmin.
{
//This opens the file the data will be recalled from; note that data_filename needs to be set
fstream file;
file.open( data_filename.str().c_str() , fstream::in);
//read first line and close file
string header_string;
getline(data_file,header_string);
stringstream header (stringstream::in | stringstream::out);
header << header_string;
file.close();
//extract data from first line
Mdouble N, t, xmin, ymin, zmin, xmax, ymax, zmax;
header >> N >> t >> xmin >> ymin >> zmin >> xmax >> ymax >> zmax;
//Set the correct formation based of dimension if the formation is not specified by the user
if (header.eof()) return 2;
else return 3;
}
| bool MD::read_next_from_data_file | ( | unsigned int | format = 0 | ) |
Reimplemented in StatisticsVector< T >.
References ParticleHandler::begin(), ParticleHandler::copyAndAddParticle(), STD_save::data_file, data_FixedParticles, dim, ParticleHandler::end(), Particle::fixParticle(), ParticleHandler::get_NumberOfParticles(), STD_save::get_options_data(), ParticleHandler::get_Particle(), STD_save::increase_counter_data(), particleHandler, ParticleHandler::removeLastParticle(), Species, t, xmax, xmin, ymax, ymin, zmax, and zmin.
Referenced by load_from_data_file(), and StatisticsVector< T >::read_next_from_data_file().
{
if (get_options_data()==2) increase_counter_data(fstream::in);
//Set the correct formation based of dimension if the formation is not specified by the user
if (format ==0)
{
switch (dim)
{
case 1:
case 2:
format=8;
break;
case 3:
format=14;
break;
}
//end case
}
//end if
static unsigned int N;
static Mdouble dummy;
//read first parameter and check if we reached the end of the file
data_file >> N;
Particle P0;
if(particleHandler.get_NumberOfParticles()<N)
while(particleHandler.get_NumberOfParticles()<N)
particleHandler.copyAndAddParticle(P0);
else
while(particleHandler.get_NumberOfParticles()>N)
particleHandler.removeLastParticle();
#ifdef DEBUG_OUTPUT
cout << "Number of particles read from file "<<N << endl;
#endif
if (data_file.eof()||data_file.peek()==-1) return false;
//read all other data available for the time step
switch(format)
{
//This is the format that Stefan's and Vitaley's code saves in - note the axis rotation
case 7:
{
data_file >> t >> xmin >> zmin >> ymin >> xmax >> zmax >> ymax;
cout << " time " << t << endl;
Mdouble radius;
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) {
data_file
>> (*it)->get_Position().X
>> (*it)->get_Position().Z
>> (*it)->get_Position().Y
>> (*it)->get_Velocity().X
>> (*it)->get_Velocity().Z
>> (*it)->get_Velocity().Y
>> radius
>> dummy;
(*it)->set_Angle(Vec3D(0.0,0.0,0.0));
(*it)->set_AngularVelocity(Vec3D(0.0,0.0,0.0));
(*it)->set_Radius(radius);
(*it)->compute_particle_mass(Species);
}
//End the interator over all particles
break;
}
//This is a 2D format
case 8:
{
data_file >> t >> xmin >> ymin >> xmax >> ymax;
zmin = 0.0;
zmax = 0.0;
Mdouble radius;
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) {
data_file
>> (*it)->get_Position().X
>> (*it)->get_Position().Y
>> (*it)->get_Velocity().X
>> (*it)->get_Velocity().Y
>> radius
>> (*it)->get_Angle().Z
>> (*it)->get_AngularVelocity().Z
>> dummy;
(*it)->get_Position().Z = 0.0;
(*it)->get_Velocity().Z = 0.0;
(*it)->get_Angle().X = 0.0;
(*it)->get_Angle().Y = 0.0;
(*it)->get_Angle().Z = -(*it)->get_Angle().Z;
(*it)->get_AngularVelocity().X = 0.0;
(*it)->get_AngularVelocity().Y = 0.0;
(*it)->get_AngularVelocity().Z = -(*it)->get_AngularVelocity().Z;
(*it)->set_Radius(radius);
(*it)->compute_particle_mass(Species);
} //end for all particles
break;
}
//This is a 3D format
case 14:
{
data_file >> t >> xmin >> ymin >> zmin >> xmax >> ymax >> zmax;
Mdouble radius;
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) {
data_file
>> (*it)->get_Position()
>> (*it)->get_Velocity()
>> radius
>> (*it)->get_Angle()
>> (*it)->get_AngularVelocity()
>> dummy;
(*it)->set_Radius(radius);
(*it)->compute_particle_mass(Species);
} //end for all particles
break;
}
//This is a 3D format
case 15:
{
data_file >> t >> xmin >> ymin >> zmin >> xmax >> ymax >> zmax;
cout<<t<<" "<<xmin<<" "<<ymin<<" "<<zmin<<" "<<xmax<<" "<<ymax<<" "<<zmax<<endl;
Mdouble radius;
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) {
data_file
>> (*it)->get_Position()
>> (*it)->get_Velocity()
>> radius
>> (*it)->get_Angle()
>> (*it)->get_AngularVelocity()
>> dummy >> dummy;
(*it)->set_Radius(radius);
(*it)->compute_particle_mass(Species);
} //end for all particles
break;
}
} //end switch statement
//fix particles, if data_FixedParticles!=0
if (data_FixedParticles) {
//cout << "fix " << min(data_FixedParticles,get_ParticleHandler().get_NumberOfParticles()) << " Particles" << endl;
for (int i=0; i<min(data_FixedParticles,(int) particleHandler.get_NumberOfParticles()); i++)
particleHandler.get_Particle(i)->fixParticle();
}
//clean up tangential springs
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it)
(*it)->get_TangentialSprings().resize(0);
return true;
}
| void MD::read_v1 | ( | std::istream & | is | ) | [virtual] |
Reads all MD data.
References Particle::compute_particle_mass(), ParticleHandler::copyAndAddParticle(), dt, gravity, max_radius, STD_save::options_data, STD_save::options_ene, STD_save::options_fstat, particleHandler, STD_save::problem_name, save_count_data, save_count_ene, save_count_fstat, save_count_stat, Species, t, tmax, Walls, WallsPeriodic, xmax, xmin, ymax, ymin, zmax, and zmin.
Referenced by read().
{
cout << "reading restart data version 1" << endl;
is >> gravity
>> xmin >> xmax >> ymin >> ymax >> zmin >> zmax
>> dt >> t >> tmax
>> save_count_data >> max_radius;
save_count_ene=save_count_data;
save_count_fstat=save_count_data;
save_count_stat=save_count_data;
string str;
is >> str;
problem_name.str(str);
int N;
unsigned int NParticles;
Particle p0;
is >> N; Species.resize(N);
is >> NParticles;
is >> N; Walls.resize(N);
is >> N; WallsPeriodic.resize(N);
is >> options_fstat >> options_data >> options_ene;
for (vector<CSpecies>::iterator it = Species.begin(); it!=Species.end(); ++it) is >> (*it);
for (unsigned int i=0;i<NParticles;i++) { is>>p0;p0.compute_particle_mass(Species);particleHandler.copyAndAddParticle(p0);}
for (vector<CWall>::iterator it = Walls.begin(); it!=Walls.end(); ++it) is >> (*it);
for (vector<CWallPeriodic>::iterator it = WallsPeriodic.begin(); it!=WallsPeriodic.end(); ++it) is >> (*it);
}
| void MD::read_v2 | ( | std::istream & | is | ) | [virtual] |
References Particle::compute_particle_mass(), ParticleHandler::copyAndAddParticle(), data_FixedParticles, dim, dt, gravity, STD_save::options_data, STD_save::options_ene, STD_save::options_fstat, STD_save::options_restart, particleHandler, save_count_data, save_count_ene, save_count_fstat, save_count_stat, Species, t, tmax, Walls, WallsPeriodic, xmax, xmin, ymax, ymin, zmax, and zmin.
Referenced by read().
{
string dummy;
cout << "version 2" << endl;
is >> dummy >> xmin
>> dummy >> xmax
>> dummy >> ymin
>> dummy >> ymax
>> dummy >> zmin
>> dummy >> zmax
>> dummy >> dt
>> dummy >> t
>> dummy >> tmax
>> dummy >> save_count_data
>> dummy >> dim
>> dummy >> gravity
>> dummy >> options_fstat
>> dummy >> options_data
>> dummy >> options_ene
>> dummy;
save_count_ene=save_count_data;
save_count_fstat=save_count_data;
save_count_stat=save_count_data;
//this is optional to allow restart files with and without options_restart
if (!strcmp(dummy.c_str(),"options_restart")) is >> options_restart >> dummy;
unsigned int N;
is >> N;
Species.resize(N);
//please don't change this code to iterators; the mixed species requires indexing
//old code: for (vector<CSpecies>::iterator it = Species.begin(); it!=Species.end(); ++it) it->read(is);
for (unsigned int i=0; i<N; i++) {
Species[i].read(is);
Species[i].MixedSpecies.resize(i);
for (unsigned int j=0; j<i; j++) {
Species[i].MixedSpecies[j].read(is);
}
}
is >> dummy >> N;
Walls.resize(N);
for (vector<CWall>::iterator it = Walls.begin(); it!=Walls.end(); ++it) it->read(is);
is >> dummy >> N;
WallsPeriodic.resize(N);
for (vector<CWallPeriodic>::iterator it = WallsPeriodic.begin(); it!=WallsPeriodic.end(); ++it) it->read(is);
is >> dummy >> N;
Particle p0;
//for each line of particle data, the line is first completely read. This is to read restart datafiles both with and without tangential spring data.
string line_string;
getline(is,line_string);
for (unsigned int i=0;i<N;i++) {
getline(is,line_string);
stringstream line (stringstream::in | stringstream::out);
line << line_string;
line >> p0;
p0.compute_particle_mass(Species);
particleHandler.copyAndAddParticle(p0);
}
//After the lines storing particle data, an optional variable 'FixedParticles' can be read in, which fixes the first 'FixedParticles' particles
if (is.peek()==70) is >> dummy >> data_FixedParticles;
}
| int MD::readArguments | ( | unsigned int | argc, |
| char * | argv[] | ||
| ) |
Can interpret main function input arguments that are passed by the driver codes.
References readNextArgument().
{
bool isRead = true;
for (unsigned int i = 1; i<argc; i+=2) {
cout << "interpreting input argument " << argv[i];
for (unsigned int j = i+1; j<argc; j++) {
if (argv[j][0]=='-') break;
cout << " " << argv[j];
}
cout << endl;
isRead &= readNextArgument(i, argc, argv);
if (!isRead) {
cerr << "Warning: not all arguments read correctly!" << endl;
exit(-10);
}
}
return isRead;
}
| int MD::readNextArgument | ( | unsigned int & | i, |
| unsigned int | argc, | ||
| char * | argv[] | ||
| ) | [virtual] |
argv[i+1] interpreted as argument of type char*, Mdouble, integer or boolean unless noted
-gravity requires three arguments
-restart or -r loads a restart file. By default, it loads <get_name()>.restart. If an argument <arg> is given it loads the file <arg>, or <arg>.restart (if the ending is not given).
Reimplemented in ChuteWithHopper, and Chute.
References get_disp(), get_dispt(), get_dt(), get_k(), get_kt(), get_mu(), STD_save::get_name(), get_tmax(), load_restart_data(), random, RNG::randomise(), set_append(), STD_save::set_counter(), set_dim(), set_dim_particle(), set_disp(), set_dispt(), set_dt(), set_FixedParticles(), set_gravity(), set_Hertzian(), set_k(), set_kt(), set_mu(), set_name(), STD_save::set_options_data(), STD_save::set_options_ene(), STD_save::set_options_fstat(), STD_save::set_options_restart(), STD_save::set_options_stat(), set_rho(), set_savecount(), set_tmax(), set_xmax(), set_xmin(), set_ymax(), set_ymin(), set_zmax(), and set_zmin().
Referenced by readArguments().
{
if (!strcmp(argv[i],"-name")) {
set_name(argv[i+1]);
} else if (!strcmp(argv[i],"-xmin")) {
set_xmin(atof(argv[i+1]));
} else if (!strcmp(argv[i],"-ymin")) {
set_ymin(atof(argv[i+1]));
} else if (!strcmp(argv[i],"-zmin")) {
set_zmin(atof(argv[i+1]));
} else if (!strcmp(argv[i],"-xmax")) {
set_xmax(atof(argv[i+1]));
} else if (!strcmp(argv[i],"-ymax")) {
set_ymax(atof(argv[i+1]));
} else if (!strcmp(argv[i],"-zmax")) {
set_zmax(atof(argv[i+1]));
} else if (!strcmp(argv[i],"-dt")) {
Mdouble old=get_dt();
set_dt(atof(argv[i+1]));
cout << " reset dt from " << old << " to " << get_dt() << endl;
} else if (!strcmp(argv[i],"-Hertzian")) {
if (i+1>=argc||argv[i+1][0]=='-') set_Hertzian(true);
else set_Hertzian(atoi(argv[i+1]));
} else if (!strcmp(argv[i],"-tmax")) {
Mdouble old = get_tmax();
set_tmax(atof(argv[i+1]));
cout << " reset tmax from " << old << " to " << get_tmax() << endl;
} else if (!strcmp(argv[i],"-save_count") || !strcmp(argv[i],"-savecount")) {
set_savecount(atoi(argv[i+1]));
} else if (!strcmp(argv[i],"-dim")) {
set_dim(atoi(argv[i+1]));
} else if (!strcmp(argv[i],"-gravity")) {
set_gravity(Vec3D(atof(argv[i+1]),atof(argv[i+2]),atof(argv[i+3])));
i+=2;
} else if (!strcmp(argv[i],"-options_fstat")) {
set_options_fstat(atoi(argv[i+1]));
} else if (!strcmp(argv[i],"-options_restart")) {
set_options_restart(atoi(argv[i+1]));
} else if (!strcmp(argv[i],"-options_data")) {
set_options_data(atoi(argv[i+1]));
} else if (!strcmp(argv[i],"-options_stat")) {
set_options_stat(atoi(argv[i+1]));
} else if (!strcmp(argv[i],"-options_ene")) {
set_options_ene(atoi(argv[i+1]));
} else if (!strcmp(argv[i],"-restart")||!strcmp(argv[i],"-r")) {
string filename;
//use default filename if no argument is given
if (i+1>=argc||argv[i+1][0]=='-') {
i--;
filename = get_name();
cout << get_name() << endl;
} else filename = argv[i+1];
//add ".restart" if necessary
const char* fileextension = (filename.c_str()+max(0,(int)filename.length()-8));
if (strcmp(fileextension,".restart"))
filename=filename+".restart";
cout << "restart from " << filename << endl;
load_restart_data(filename);
} else if (!strcmp(argv[i],"-k")) {
Mdouble old = get_k();
set_k(atof(argv[i+1]));
cout << " reset k from " << old << " to " << get_k() << endl;
} else if (!strcmp(argv[i],"-disp")) {
Mdouble old = get_disp();
set_disp(atof(argv[i+1]));
cout << " reset disp from " << old << " to " << get_disp() << endl;
} else if (!strcmp(argv[i],"-kt")) {
Mdouble old = get_kt();
set_kt(atof(argv[i+1]));
cout << " reset kt from " << old << " to " << get_kt() << endl;
} else if (!strcmp(argv[i],"-dispt")) {
Mdouble old = get_dispt();
set_dispt(atof(argv[i+1]));
cout << " reset dispt from " << old << " to " << get_dispt() << endl;
} else if (!strcmp(argv[i],"-mu")) {
Mdouble old = get_mu();
set_mu(atof(argv[i+1]));
cout << " reset mu from " << old << " to " << get_mu() << endl;
} else if (!strcmp(argv[i],"-randomise") || !strcmp(argv[i],"-randomize")) {
random.randomise();
i--;
} else if (!strcmp(argv[i],"-append")) {
set_append(true);
i--;
} else if (!strcmp(argv[i],"-fixedParticles")) {
set_FixedParticles(atof(argv[i+1]));
} else if (!strcmp(argv[i],"-rho")) {
set_rho(atof(argv[i+1]));
} else if (!strcmp(argv[i],"-dim_particle")) {
set_dim_particle(atoi(argv[i+1]));
} else if (!strcmp(argv[i],"-counter")) {
set_counter(atoi(argv[i+1]));
} else return false;
return true; //returns true if argv is found
}
| void MD::Remove_Duplicate_Periodic_Particles | ( | ) | [private] |
Remove periodic duplicate Particles (i.e. removes particles created by Check_and_Duplicate_Periodic_Particle(int i, int nWallPeriodic)). Note that between these two functions it is not allowed to create additional functions.
Remove periodic duplicate Particles (i.e. removes particles created by Check_and_Duplicate_Periodic_Particle(int i, int nWallPeriodic)). Note that between these two functions it is not allowed to create additional functions. It returns the number of Particles removed.
References ParticleHandler::get_NumberOfParticles(), ParticleHandler::get_Particle(), Particle::get_PeriodicFromParticle(), particleHandler, PeriodicCreated, and removeParticle().
Referenced by solve(), and statistics_from_restart_data().
{
int R=0;
for(int i=particleHandler.get_NumberOfParticles()-1; i>=0&&particleHandler.get_Particle(i)->get_PeriodicFromParticle()!=NULL; i--)
{
removeParticle(i);
R++;
}
if(R!=PeriodicCreated)
{
cerr<<"ERROR :: Periodic Particles Removed not equal to Periodic Particles Created"<<endl<<endl;
exit(-1);
}
}
| void MD::Remove_Particle | ( | int | IP | ) |
This function removes partice IP from the vector of particles by moving the last particle in the vector to the position if IP Also it checks if the moved Particle has any tangentialsspring-information, which needs to be moved to a different particle, because tangential spring information always needs to be stored in the real particle with highest particle index.
| virtual void MD::removeParticle | ( | int | iP | ) | [inline, virtual] |
Referenced by Chute::clean_chute(), and Remove_Duplicate_Periodic_Particles().
| void MD::reset_DeltaMax | ( | ) | [inline, protected] |
sets the history parameter DeltaMax of all particles to zero
{for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) (*it)->get_DeltaMaxs().resize(0);}
| void MD::reset_TangentialSprings | ( | ) | [inline, protected] |
sets the history parameter TangentialSprings of all particles to zero
{for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) (*it)->get_TangentialSprings().resize(0);}
| void MD::save_restart_data | ( | ) | [virtual] |
Stores all MD data.
This function stores all MD data - See also MD::load_restart_data.
todo{This can probably done a lot nicer}
References STD_save::options_restart, STD_save::problem_name, save_restart_data_counter, and write().
Referenced by solve().
{
stringstream restart_filename;
restart_filename.str("");
restart_filename << problem_name.str().c_str() << ".restart";
if (options_restart==2) {
restart_filename << ".";
if (save_restart_data_counter<1000) restart_filename << "0";
if (save_restart_data_counter<100) restart_filename << "0";
if (save_restart_data_counter<10) restart_filename << "0";
restart_filename << save_restart_data_counter;
save_restart_data_counter++;
}
ofstream restart_data(restart_filename.str().c_str());
restart_data.precision(8); //max. 16 for Mdouble precision
if( restart_data.good() )
{
write(restart_data);
restart_data.close();
} else {
cerr << "restart_data " << restart_filename.str() << " could not be saved" << endl;
exit(-1);
}
}
| void MD::set_append | ( | bool | new_ | ) | [inline] |
| void MD::set_collision_time_and_normal_and_tangential_restitution_coefficient | ( | Mdouble | tc, |
| Mdouble | eps, | ||
| Mdouble | beta, | ||
| Mdouble | mass1, | ||
| Mdouble | mass2, | ||
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_collision_time_and_normal_and_tangential_restitution_coefficient_nodispt | ( | Mdouble | tc, |
| Mdouble | eps, | ||
| Mdouble | beta, | ||
| Mdouble | mass1, | ||
| Mdouble | mass2, | ||
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_collision_time_and_restitution_coefficient | ( | Mdouble | tc, |
| Mdouble | eps, | ||
| Mdouble | mass, | ||
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Sets k, disp such that it matches a given tc and eps for a collision of two copies of P.
| void MD::set_collision_time_and_restitution_coefficient | ( | Mdouble | tc, |
| Mdouble | eps, | ||
| Mdouble | mass1, | ||
| Mdouble | mass2, | ||
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Set k, disp such that is matches a given tc and eps for a collision of two different masses.
Recall the resitution constant is a function of k, disp and the mass of each particle in the collision See also set_collision_time_and_restitution_coefficient(Mdouble tc, Mdouble eps, Mdouble mass)
| void MD::set_depth | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_dim | ( | int | new_dim | ) | [inline] |
Allows the dimension of the simulation to be changed.
Referenced by HGRID_3D::constructor(), and readNextArgument().
{ if (new_dim>=1 && new_dim<=3) dim = new_dim; else { cerr << "Error in set_dim" << endl; exit(-1); }}
| void MD::set_dim_particle | ( | int | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Allows the dimension of the particle (f.e. for mass) to be changed.
Referenced by constructor(), HGRID_3D::constructor(), load_par_ini_file(), and readNextArgument().
| void MD::set_disp | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Allows the normal dissipation to be changed.
Referenced by load_par_ini_file(), and readNextArgument().
| void MD::set_disprolling | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_dispt | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Allows the tangential viscosity to be changed.
Referenced by readNextArgument().
| void MD::set_disptorsion | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_dissipation | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Allows the normal dissipation to be changed.
Referenced by Chute::set_collision_time_and_restitution_coefficient().
| void MD::set_do_stat_always | ( | bool | new_ | ) | [inline] |
Sets how often the data is saved using the number of saves wanted, tmax, and dt. See also set_savecount.
Referenced by constructor().
{do_stat_always=new_;}
| void MD::set_dt | ( | Mdouble | new_dt | ) | [inline] |
| void MD::set_dt | ( | Particle & | P | ) | [inline] |
Sets dt to 1/50-th of the collision time for two copies of P.
References Particle::get_IndSpecies(), and Particle::get_Mass().
{
Mdouble mass = P.get_Mass();
dt = .02 * get_collision_time(mass);
if (Species[P.get_IndSpecies()].dispt) { dt = min(dt,0.125*mass/Species[P.get_IndSpecies()].dispt); cerr << "Warning: dispt is large, dt had to be lowered"; }
}
| void MD::set_dt | ( | ) | [inline] |
Sets dt to 1/50-th of the smallest possible collision time.
Reimplemented in Chute.
Referenced by load_par_ini_file(), and readNextArgument().
{
setup_particles_initial_conditions();
for (unsigned int i=0;i<particleHandler.get_NumberOfParticles();i++) particleHandler.get_Particle(i)->compute_particle_mass(Species);
set_dt(get_SmallestParticle()[0]);
}
| void MD::set_dt_by_mass | ( | Mdouble | mass | ) | [inline] |
Sets dt to 1/50-th of the collision time for two particles of mass P.
{dt = .02 * get_collision_time(mass);}
| void MD::set_ene_ela | ( | Mdouble | new_ | ) | [inline] |
Sets ene_ela.
{ene_ela=new_;}
| void MD::set_FixedParticles | ( | unsigned int | n | ) | [inline, protected] |
Referenced by readNextArgument().
{for (unsigned int i=0; i<std::min((unsigned int)particleHandler.get_NumberOfParticles(),n); i++) particleHandler.get_Particle(i)->fixParticle();}
| void MD::set_format | ( | int | new_ | ) | [inline] |
{format = new_;};
| void MD::set_gravity | ( | Vec3D | new_gravity | ) | [inline] |
Allows the gravitational acceleration to be changed.
Referenced by load_par_ini_file(), readNextArgument(), and Chute::set_ChuteAngle().
{gravity = new_gravity;}
| void MD::set_Hertzian | ( | bool | new_ | ) | [inline] |
Referenced by readNextArgument().
| virtual void MD::set_initial_pressures_for_pressure_controlled_walls | ( | ) | [inline, protected, virtual] |
Referenced by solve().
{};
Allows the spring constant to be changed.
Referenced by constructor(), load_par_ini_file(), readNextArgument(), and Chute::set_collision_time_and_restitution_coefficient().
| void MD::set_k1 | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_k2max | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_k_and_restitution_coefficient | ( | Mdouble | k_, |
| Mdouble | eps, | ||
| Mdouble | mass, | ||
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_kc | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_krolling | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_kt | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Allows the spring constant to be changed.
Referenced by readNextArgument().
| void MD::set_ktorsion | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_MixedSpecies | ( | int | i, |
| int | j, | ||
| CSpecies & | S | ||
| ) | [inline] |
Allows the mixed species to be set.
References CSpecies::MixedSpecies.
| void MD::set_mu | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Allows the Coulomb friction coefficient to be changed.
Referenced by ChuteBottom::make_rough_bottom(), and readNextArgument().
| void MD::set_murolling | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_mutorsion | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_name | ( | const char * | name | ) | [inline] |
Sets the name of the problem, used for the same data files.
Reimplemented from STD_save.
Referenced by ChuteBottom::constructor(), and readNextArgument().
{problem_name.str(""); problem_name << name;}
| void MD::set_number_of_saves | ( | Mdouble | N | ) | [inline] |
Referenced by Chute::readNextArgument().
{set_number_of_saves_all(N);}
| void MD::set_number_of_saves_all | ( | Mdouble | N | ) | [inline] |
{if (dt) {set_save_count_all((N>1)?((int)ceil(tmax/dt/(N-1.0))):1000000);} else {cerr << "Error in set_number_of_saves_all ; dt must be set"; exit(-1);} }
| void MD::set_number_of_saves_data | ( | Mdouble | N | ) | [inline] |
{if (dt) {save_count_data = (N>1)?((int)ceil(tmax/dt/(N-1.0))):1000000; } else {cerr << "Error in set_number_of_saves_data ; dt must be set"; exit(-1);} }
| void MD::set_number_of_saves_ene | ( | Mdouble | N | ) | [inline] |
{if (dt) {save_count_ene = (N>1)?((int)ceil(tmax/dt/(N-1.0))):1000000; } else {cerr << "Error in set_number_of_saves_ene ; dt must be set"; exit(-1);} }
| void MD::set_number_of_saves_fstat | ( | Mdouble | N | ) | [inline] |
{if (dt) {save_count_fstat = (N>1)?((int)ceil(tmax/dt/(N-1.0))):1000000; } else {cerr << "Error in set_number_of_saves_fstat; dt must be set"; exit(-1);} }
| void MD::set_number_of_saves_stat | ( | Mdouble | N | ) | [inline] |
{if (dt) {save_count_stat = (N>1)?((int)ceil(tmax/dt/(N-1.0))):1000000; } else {cerr << "Error in set_number_of_saves_stat ; dt must be set"; exit(-1);} }
| void MD::set_NWall | ( | int | N | ) | [inline] |
Allows the number of walls to be changed.
Referenced by ChuteWithHopper::add_hopper(), constructor(), load_par_ini_file(), and Chute::setup_particles_initial_conditions().
{if (N>=0) {Walls.resize(N);} else {cerr << "Error in set_NWall" << endl; exit(-1);}}
| void MD::set_NWallPeriodic | ( | int | N | ) | [inline] |
Allows the number of periodic walls to be changed.
Referenced by constructor(), load_par_ini_file(), and Chute::setup_particles_initial_conditions().
{if (N>=0) {WallsPeriodic.resize(N);} else {cerr << "Error in set_NWallPeriodic" << endl; exit(-1);}}
| void MD::set_plastic_k1_k2max_kc_depth | ( | Mdouble | k1_, |
| Mdouble | k2max_, | ||
| Mdouble | kc_, | ||
| Mdouble | depth_, | ||
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
| void MD::set_restart_version | ( | int | new_ | ) | [inline] |
Sets restart_version.
{restart_version=new_;}
| void MD::set_restarted | ( | bool | new_ | ) | [inline] |
Referenced by constructor(), and load_restart_data().
{
restarted=new_;
set_append(new_);
}
| void MD::set_rho | ( | Mdouble | new_, |
| unsigned int | indSpecies = 0 |
||
| ) | [inline] |
Allows the density to be changed.
Referenced by constructor(), load_par_ini_file(), and readNextArgument().
| void MD::set_save_count_all | ( | int | new_ | ) | [inline] |
Referenced by constructor().
{if (new_>0) {save_count_data =new_;save_count_ene =new_;save_count_stat =new_;save_count_fstat=new_;} else {cerr << "Error in set_save_count_all (set_save_count_all("<<new_<<"))"<< endl; exit(-1);}}
| void MD::set_save_count_data | ( | int | new_ | ) | [inline] |
Referenced by load_par_ini_file().
{if (new_>0) {save_count_data =new_;} else {cerr << "Error in set_save_count_data, (set_save_count_data ("<<new_<<"))"<< endl; exit(-1);}}
| void MD::set_save_count_ene | ( | int | new_ | ) | [inline] |
{if (new_>0) {save_count_ene =new_;} else {cerr << "Error in set_save_count_ene, (set_save_count_ene ("<<new_<<"))"<< endl; exit(-1);}}
| void MD::set_save_count_fstat | ( | int | new_ | ) | [inline] |
Referenced by load_par_ini_file().
{if (new_>0) {save_count_fstat=new_;} else {cerr << "Error in set_save_count_fstat, (set_save_count_fstat("<<new_<<"))"<< endl; exit(-1);}}
| void MD::set_save_count_stat | ( | int | new_ | ) | [inline] |
{if (new_>0) {save_count_stat =new_;} else {cerr << "Error in set_save_count_stat, (set_save_count_stat ("<<new_<<"))"<< endl; exit(-1);}}
| void MD::set_savecount | ( | int | new_ | ) | [inline] |
Allows the number of time steps between saves to be changed, see also set_number_of_saves.
Referenced by load_par_ini_file(), ChuteBottom::make_rough_bottom(), and readNextArgument().
{if (new_>0) {set_save_count_all(new_);} else {cerr << "Error in set_savecount (set_savecount("<<new_<<"))"<< endl; exit(-1);}}
| void MD::set_tmax | ( | Mdouble | new_tmax | ) | [inline] |
Allows the upper time limit to be changed.
Referenced by load_par_ini_file(), ChuteBottom::make_rough_bottom(), and readNextArgument().
{if (new_tmax>=0){tmax=new_tmax;} else { cerr << "Error in set_tmax" << endl; exit(-1); }}
| void MD::set_xballs_additional_arguments | ( | string | new_ | ) | [inline] |
Set the additional arguments for xballs.
{xballs_additional_arguments=new_.c_str();}
| void MD::set_xballs_cmode | ( | int | new_cmode | ) | [inline] |
{xballs_cmode=new_cmode;}
| void MD::set_xballs_colour_mode | ( | int | new_cmode | ) | [inline] |
Set the xball output mode.
{xballs_cmode=new_cmode;}
| void MD::set_xballs_scale | ( | Mdouble | new_scale | ) | [inline] |
Set the scale of the xballs problem. The default is fit to screen.
{xballs_scale=new_scale;}
| void MD::set_xballs_vector_scale | ( | double | new_vscale | ) | [inline] |
Set the scale of vectors in xballs.
{xballs_vscale=new_vscale;}
| void MD::set_xmax | ( | Mdouble | new_xmax | ) | [inline] |
Sets xmax and walls, assuming the standard definition of walls as in the default constructor.
Referenced by readNextArgument(), Chute::set_ChuteLength(), ChuteWithHopper::set_ChuteLength(), and ChuteWithHopper::set_shift().
| void MD::set_xmin | ( | Mdouble | new_xmin | ) | [inline] |
Sets xmin and walls, assuming the standard definition of walls as in the default constructor.
Referenced by readNextArgument(), and ChuteWithHopper::set_ChuteLength().
| void MD::set_ymax | ( | Mdouble | new_ymax | ) | [inline] |
Sets ymax and walls, assuming the standard definition of walls as in the default constructor.
Referenced by readNextArgument(), and Chute::set_ChuteWidth().
| void MD::set_ymin | ( | Mdouble | new_ymin | ) | [inline] |
Referenced by readNextArgument().
| void MD::set_zmax | ( | Mdouble | new_zmax | ) | [inline] |
Sets ymax and walls, assuming the standard definition of walls as in the default constructor.
Referenced by ChuteWithHopper::add_hopper(), Chute::readNextArgument(), readNextArgument(), and Chute::set_InflowHeight().
| void MD::set_zmin | ( | Mdouble | new_zmin | ) | [inline] |
Sets ymin and walls, assuming the standard definition of walls as in the default constructor.
Referenced by readNextArgument().
| void MD::setup_particles_initial_conditions | ( | ) | [virtual] |
This function allows the initial conditions of the particles to be set, by default locations is random.
This setup the particles initial conditions it is virtual as you expect the user to override this.
Remember particle properties must also be defined here.
By default the particles are randomly disibuted///
Reimplemented in ChuteBottom, ChuteWithHopper, and Chute.
Referenced by solve().
{
/*Particle p0;
for (unsigned int i=0;i<2;i++)
{
p0.set_Position(Vec3D(random.get_RN(0.0,get_xmax()),random.get_RN(0.0,get_ymax()),0.0));
p0.set_Velocity(Vec3D(0.0,0.0,0.0));
p0.set_Radius(0.0005);
p0.compute_particle_mass(Species);
particleHandler.copyAndAddParticle(p0);
} //end loop over partilces*/
}
| void MD::solve | ( | ) |
The work horse of the code.
This is the main solve loop where all the action takes place///.
This is where the main code works
Print counter if it is used
Set up the data file names
Initialise the time and sets up the initial conditions for the simulation
Set up the data file names
This creates the file the data will be saved to
Set up the ene file names
This creates the file ene will be saved to
Set up the fstat file names
This creates the file fstatistics will be saved to
initializes PrevPosition if PositionVerlet is used
Setup the mass of each particle.
Other initializations
velocity verlet requires initial force calculation
This is the main loop over advancing time
Here we output to the data file
Check if rotation is on
Loop over all particles doing the time integration step
Output statistical information
Compute forces
Loop over all particles doing the time integration step
References actions_after_solve(), actions_after_time_step(), actions_before_time_loop(), actions_before_time_step(), append, ParticleHandler::begin(), Check_and_Duplicate_Periodic_Particles(), compute_all_forces(), compute_particle_masses(), continue_solve(), cout_time(), create_xballs_script(), STD_save::data_file, do_integration_after_force_computation(), do_integration_before_force_computation(), do_stat_always, dt, ParticleHandler::end(), ene_ela, STD_save::ene_file, finish_statistics(), STD_save::fstat_file, fstat_header(), get_append(), STD_save::get_counter(), get_LargestParticle(), get_mu(), STD_save::get_options_data(), STD_save::get_options_ene(), STD_save::get_options_fstat(), STD_save::get_options_restart(), STD_save::get_options_stat(), ParticleHandler::get_Particle(), get_ParticleHandler(), Particle::get_Position(), Particle::get_Radius(), Particle::get_Velocity(), gravity, HGRID_actions_after_integration(), HGRID_actions_before_integration(), HGRID_actions_before_time_loop(), HGRID_actions_before_time_step(), STD_save::increase_counter_data(), STD_save::increase_counter_ene(), STD_save::increase_counter_fstat(), STD_save::increase_counter_stat(), initialize_statistics(), max_radius, STD_save::open_data_file(), STD_save::open_ene_file(), STD_save::open_fstat_file(), output_ene(), output_statistics(), output_xballs_data(), particleHandler, PeriodicCreated, STD_save::problem_name, process_statistics(), Remove_Duplicate_Periodic_Particles(), restarted, rotation, save_count_data, save_count_ene, save_count_fstat, save_count_stat, save_data_data, save_data_ene, save_data_fstat, save_data_stat, save_restart_data(), STD_save::set_data_filename(), STD_save::set_ene_filename(), STD_save::set_fstat_filename(), set_initial_pressures_for_pressure_controlled_walls(), setup_particles_initial_conditions(), Species, start_ene(), t, and tmax.
Referenced by ChuteBottom::make_rough_bottom().
{
#ifdef DEBUG_OUTPUT
cerr << "Entered solve" << endl;
#endif
if (get_counter()) cout << "Count: " << get_counter() << endl;
else cout << "Started to solve ..." << endl;
if (get_counter()>0) problem_name << "." << get_counter();
if (!restarted) {
t=0;
setup_particles_initial_conditions();
#ifdef DEBUG_OUTPUT
cerr << "Have created the particles initial conditions " << endl;
#endif
}
set_data_filename();
//undo append if data file does not exist
//~ cout << " get_data_filename()" << get_data_filename() << ", append" << append << ", access" << access(get_data_filename().c_str(), F_OK ) << endl;
//~ if (append && (access(get_data_filename().c_str(), F_OK ) == -1) ) {
//~ append = false;
//~ cerr << "Warning: cannot append because data file does not exist; set append to false." << endl;
//~ }
if (get_append()) {
if (get_options_data()==2) increase_counter_data(fstream::out);
if (get_options_ene()==2) increase_counter_ene(fstream::out);
if (get_options_stat()==2) increase_counter_stat(fstream::out);
if (get_options_fstat()==2) increase_counter_fstat(fstream::out);
}
fstream::openmode mode;
if (append) mode=fstream::out|fstream::app;
else mode=fstream::out;
if(!open_data_file(mode)) {cerr << "Problem opening data file aborting" << endl; exit(-1);}
set_ene_filename();
if(!open_ene_file(mode)){cerr << "Problem opening ene file aborting" << endl; exit(-1);}
set_fstat_filename();
if(!open_fstat_file(mode)){cerr << "Problem opening fstat file aborting" << endl; exit(-1);}
#ifdef USE_SIMPLE_VERLET_INTEGRATION
cout << "using simple verlet integration" << endl;
for (int i = 0;i<Particles.size();i++)
get_ParticleHandler().get_Particle(i)->PrevPosition = get_ParticleHandler().get_Particle(i)->get_Position() - get_ParticleHandler().get_Particle(i)->get_Velocity() * dt + 0.5 * gravity * dt * dt;
#endif
compute_particle_masses();
max_radius = get_LargestParticle()->get_Radius();
actions_before_time_loop();
initialize_statistics();
HGRID_actions_before_time_loop();
ene_ela = 0;
if (!append) start_ene();
//needed for output_statistics
save_data_data =false;
save_data_ene =false;
save_data_stat =false;
save_data_fstat=false;
if(do_stat_always||save_data_stat) output_statistics();
PeriodicCreated=Check_and_Duplicate_Periodic_Particles();
HGRID_actions_before_time_step();
Mdouble dt_=dt; dt=0; compute_all_forces(); dt=dt_;
Remove_Duplicate_Periodic_Particles();
set_initial_pressures_for_pressure_controlled_walls();
if(do_stat_always||save_data_stat) process_statistics(save_data_stat);
#ifdef DEBUG_OUTPUT
cerr << "Have computed the initial values for the forces " << endl;
#endif
if (append) {
save_data_data =false;
save_data_ene =false;
save_data_stat =false;
save_data_fstat=false;
} else {
save_data_data =true;
save_data_ene =true;
save_data_stat =true;
save_data_fstat=true;
}
int count_data=1;
int count_ene=1;
int count_stat=1;
int count_fstat=1;
//create .disp file if .data file is created
if (get_options_data()) create_xballs_script();
ene_ela = 0;
while (t-dt<tmax&&continue_solve())
{
if (save_data_data)
{
if (get_options_data()==2) increase_counter_data(fstream::out);
if (get_options_restart()) save_restart_data();
}
if (save_data_ene)
{
if (get_options_ene()==2) increase_counter_ene(fstream::out);
}
if(save_data_stat)
{
if (get_options_stat()==2) increase_counter_stat(fstream::out);
}
if (save_data_fstat)
{
if (get_options_fstat()==2) increase_counter_fstat(fstream::out);
}
rotation=false;
for (unsigned int i=0;i<Species.size();i++) {
if (Species[i].get_mu()||Species[i].get_murolling()||Species[i].get_mutorsion())
rotation = true;
for (unsigned int j=0;j<Species[i].MixedSpecies.size();j++) {
if (Species[i].MixedSpecies[j].get_mu()||Species[i].MixedSpecies[j].get_murolling()||Species[i].MixedSpecies[j].get_mutorsion())
rotation = true;
}
}
HGRID_actions_before_integration();
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it)
{
do_integration_before_force_computation(*it);
}
HGRID_actions_after_integration();
if (save_data_data) output_xballs_data();
if (save_data_stat||do_stat_always) output_statistics();
actions_before_time_step();
PeriodicCreated=Check_and_Duplicate_Periodic_Particles();
HGRID_actions_before_time_step();
if (save_data_fstat) fstat_header();
compute_all_forces();
Remove_Duplicate_Periodic_Particles();
actions_after_time_step();
HGRID_actions_before_integration();
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it)
{
do_integration_after_force_computation(*it);
}
HGRID_actions_after_integration();
if (save_data_ene) output_ene();
if(save_data_stat||do_stat_always) process_statistics(save_data_stat);
t+=dt;
//To write the last time step (acutaly data is written at t+0.5*dt where t alreay is larger then tmax)
if (count_data == save_count_data || t>tmax) //write data
{
save_data_data=true;
count_data=1;
cout_time();
}
else
{
save_data_data=false;
count_data++;
}
if (count_ene == save_count_ene || t>tmax) //write data
{
save_data_ene=true;
count_ene=1;
}
else
{
save_data_ene=false;
count_ene++;
}
if (count_fstat == save_count_fstat || t>tmax) //write data
{
save_data_fstat=true;
count_fstat=1;
}
else
{
save_data_fstat=false;
count_fstat++;
}
if (count_stat == save_count_stat || t>tmax) //write data
{
save_data_stat=true;
count_stat=1;
}
else
{
save_data_stat=false;
count_stat++;
}
}
//end loop over interaction count
actions_after_solve();
cout << endl;
//To make sure get_t gets the correct time for outputting statistics
finish_statistics();
data_file.close();
ene_file.close();
fstat_file.close();
}
| void MD::solve | ( | unsigned int | argc, |
| char * | argv[] | ||
| ) | [inline] |
Read arguments before solving.
{
readArguments(argc, argv);
solve();
};
| void MD::solveWithMDCLR | ( | ) |
Tries to solve the problem using MDCLR.
| void MD::start_ene | ( | ) | [protected, virtual] |
Functions for ene file.
This function gathers statistical data.
todo{Why is there a +6 here?}
References append, and STD_save::ene_file.
Referenced by solve(), and statistics_from_restart_data().
{
//only write if we don't restart
if (append) return;
static int width = ene_file.precision() + 6;
ene_file
<< setw(width) << "t" << " "
<< setw(width) << "ene_gra" << " "
<< setw(width) << "ene_kin" << " "
<< setw(width) << "ene_rot" << " "
<< setw(width) << "ene_ela" << " "
<< setw(width) << "X_COM" << " "
<< setw(width) << "Y_COM" << " "
<< setw(width) << "Z_COM" << endl;
}
| void MD::statistics_from_restart_data | ( | const char * | name | ) |
Loads all MD data and plots statistics for all timesteps in the .data file.
statistics_from_restart_data
todo{Check this whole function}
References actions_after_time_step(), actions_before_time_loop(), actions_before_time_step(), Check_and_Duplicate_Periodic_Particles(), compute_all_forces(), compute_particle_masses(), STD_save::data_file, STD_save::data_filename, HGRID_actions_before_time_loop(), HGRID_actions_before_time_step(), initialize_statistics(), load_from_data_file(), load_restart_data(), output_ene(), PeriodicCreated, Remove_Duplicate_Periodic_Particles(), save_data_data, save_data_ene, save_data_fstat, save_data_stat, start_ene(), STD_save::stat_file, STD_save::stat_filename, and t.
{
//This function loads all MD data
load_restart_data();
//This creates the file statistics will be saved to
stat_filename.str("");
stat_filename << name << ".stat";
stat_file.open( stat_filename.str().c_str() , fstream::out);
// Sets up the initial conditions for the simulation
// setup_particles_initial_conditions();
// Setup the previous position arrays and mass of each particle.
compute_particle_masses();
// Other routines required to jump-start the simulation
actions_before_time_loop();
initialize_statistics();
HGRID_actions_before_time_loop();
start_ene();
while (load_from_data_file(data_filename.str().c_str())) {
HGRID_actions_before_time_loop();
actions_before_time_step();
PeriodicCreated=Check_and_Duplicate_Periodic_Particles();
HGRID_actions_before_time_step();
save_data_data = true;
save_data_ene = true;
save_data_stat = true;
save_data_fstat= true;
compute_all_forces();
Remove_Duplicate_Periodic_Particles();
actions_after_time_step();
output_ene();
cout << setprecision(6) << t << endl;
}
data_file.close();
stat_file.close();
}
| void MD::write | ( | std::ostream & | os | ) | [virtual] |
Writes all MD data.
todo{the following line breaks the selftest, but should be turned on to restart properly}
Reimplemented in ChuteWithHopper, HGRID_base, and Chute.
References ParticleHandler::begin(), dim, dt, ParticleHandler::end(), ParticleHandler::get_NumberOfParticles(), gravity, STD_save::options_data, STD_save::options_ene, STD_save::options_fstat, STD_save::options_restart, particleHandler, STD_save::problem_name, save_count_data, save_count_ene, save_count_fstat, Species, t, tmax, Walls, WallsPeriodic, xmax, xmin, ymax, ymin, zmax, and zmin.
Referenced by save_restart_data().
{
os
<< "restart_version " << 3
<< endl
<< "name " << problem_name.str()
<< endl
<< "xmin " << xmin << " xmax " << xmax
<< " ymin " << ymin << " ymax " << ymax
<< " zmin " << zmin << " zmax " << zmax
<< endl
<< "dt " << dt
<< " t " << t
<< " tmax " << tmax
<< " save_count_data " << save_count_data
<< " save_count_ene " << save_count_ene
<< " save_count_stat " << save_count_fstat
<< " save_count_fstat " << save_count_fstat
<< endl
<< "dim " << dim
<< " gravity " << gravity
<< endl
<< "options_fstat " << options_fstat
<< " options_data " << options_data
<< " options_ene " << options_ene
<< " options_restart " << options_restart
<< endl;
os<< "Species " << Species.size() << endl;
for (vector<CSpecies>::iterator it = Species.begin(); it!=Species.end(); ++it) {
os << (*it) << endl;
for (vector<CSpecies>::iterator it2 = it->MixedSpecies.begin(); it2!=it->MixedSpecies.end(); ++it2) os << (*it2) << " (mixed)" << endl;
}
os<< "Walls " << Walls.size() << endl;
for (vector<CWall>::iterator it = Walls.begin(); it!=Walls.end(); ++it) os << (*it) << endl;
os<< "WallsPeriodic " << WallsPeriodic.size() << endl;
for (vector<CWallPeriodic>::iterator it = WallsPeriodic.begin(); it!=WallsPeriodic.end(); ++it) os << (*it) << endl;
os<< "Particles " << particleHandler.get_NumberOfParticles() << endl;
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) {
os << (**it);
if (Species.size()>1) os << (*it)->get_IndSpecies();
os << endl;
}
//should be optional: os << "FixedParticles " << data_FixedParticles << endl;
}
| void MD::write_v1 | ( | std::ostream & | os | ) | [virtual] |
Writes all MD data.
References ParticleHandler::begin(), dim, dt, ParticleHandler::end(), ParticleHandler::get_NumberOfParticles(), gravity, max_radius, STD_save::options_data, STD_save::options_ene, STD_save::options_fstat, particleHandler, STD_save::problem_name, save_count_data, Species, t, tmax, Walls, WallsPeriodic, xmax, xmin, ymax, ymin, zmax, and zmin.
{
os << dim << " " << gravity << " "
<< xmin << " " << xmax << " " << ymin << " " << ymax << " " << zmin << " " << zmax << " "
<< dt << " " << t << " " << tmax << " "
<< save_count_data << " " << max_radius << " "
<< problem_name.str() << " "
<< Species.size() << " "
<< particleHandler.get_NumberOfParticles() << " "
<< Walls.size() << " "
<< WallsPeriodic.size() << endl;
os << options_fstat << " " << options_data << " " << options_ene << endl;
for (vector<CSpecies>::iterator it = Species.begin(); it!=Species.end(); ++it) os << (*it) << endl;
for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) os << (*it) << endl;
for (vector<CWall>::iterator it = Walls.begin(); it!=Walls.end(); ++it) os << (*it) << endl;
for (vector<CWallPeriodic>::iterator it = WallsPeriodic.begin(); it!=WallsPeriodic.end(); ++it) os << (*it) << endl;
}
bool MD::append [private] |
Referenced by solve(), and start_ene().
int MD::data_FixedParticles [private] |
Referenced by constructor(), read(), read_next_from_data_file(), and read_v2().
int MD::dim [private] |
The dimension of the simulation.
Referenced by constructor(), output_xballs_data(), print(), read(), read_next_from_data_file(), read_v2(), write(), and write_v1().
bool MD::do_stat_always [private] |
Referenced by compute_internal_forces(), compute_plastic_internal_forces(), compute_walls(), and solve().
These are numerical constants like the time step size.
Referenced by compute_internal_forces(), compute_plastic_internal_forces(), compute_walls(), constructor(), do_integration_after_force_computation(), do_integration_before_force_computation(), print(), read(), read_v1(), read_v2(), solve(), write(), and write_v1().
Mdouble MD::ene_ela [private] |
used in force calculations
Referenced by compute_internal_forces(), compute_plastic_internal_forces(), compute_walls(), output_ene(), and solve().
int MD::format [private] |
Reimplemented in StatisticsVector< T >.
Referenced by constructor(), and output_xballs_data().
Vec3D MD::gravity [private] |
Gravitational acceleration.
Referenced by compute_external_forces(), constructor(), output_ene(), print(), read(), read_v1(), read_v2(), Chute::set_ChuteAngle(), solve(), write(), and write_v1().
bool MD::Hertzian [private] |
Referenced by compute_plastic_internal_forces(), and constructor().
Mdouble MD::max_radius [private] |
Referenced by Check_and_Duplicate_Periodic_Particle(), read_v1(), solve(), and write_v1().
ParticleHandler MD::particleHandler [private] |
Referenced by Check_and_Duplicate_Periodic_Particle(), Check_and_Duplicate_Periodic_Particles(), compute_all_forces(), constructor(), output_ene(), output_xballs_data(), print(), read(), read_next_from_data_file(), read_v1(), read_v2(), Remove_Duplicate_Periodic_Particles(), solve(), write(), and write_v1().
int MD::PeriodicCreated [private] |
Referenced by constructor(), Remove_Duplicate_Periodic_Particles(), solve(), and statistics_from_restart_data().
int MD::restart_version [private] |
Referenced by read().
bool MD::restarted [private] |
Referenced by solve().
bool MD::rotation [private] |
Referenced by do_integration_after_force_computation(), do_integration_before_force_computation(), and solve().
int MD::save_count_data [private] |
int MD::save_count_ene [private] |
int MD::save_count_fstat [private] |
bool MD::save_data_data [private] |
Referenced by solve(), and statistics_from_restart_data().
bool MD::save_data_ene [private] |
bool MD::save_data_fstat [private] |
bool MD::save_data_stat [private] |
int MD::save_restart_data_counter [private] |
Referenced by constructor(), and save_restart_data().
vector<CSpecies> MD::Species [protected] |
These are the particle parameters like dissipation etc.
Referenced by add_Species(), compute_internal_forces(), compute_plastic_internal_forces(), compute_walls(), constructor(), Chute::create_inflow_particle(), ChuteWithHopper::create_inflow_particle(), Chute::get_collision_time(), Chute::get_LargestParticle(), Chute::get_restitution_coefficient(), Chute::get_SmallestParticle(), Chute::initialize_inflow_particle(), print(), read(), read_next_from_data_file(), read_v1(), read_v2(), Chute::set_collision_time_and_restitution_coefficient(), ChuteBottom::setup_particles_initial_conditions(), solve(), write(), and write_v1().
This stores the current time.
Referenced by compute_internal_forces(), compute_plastic_internal_forces(), compute_walls(), find_next_data_file(), output_xballs_data(), print(), read(), read_dim_from_data_file(), read_next_from_data_file(), read_v1(), read_v2(), StatisticsVector< T >::set_tintStat(), StatisticsVector< T >::set_tmaxStat(), StatisticsVector< T >::set_tminStat(), solve(), statistics_from_restart_data(), write(), and write_v1().
Referenced by constructor(), print(), read(), read_v1(), read_v2(), solve(), write(), and write_v1().
This is a particle class, which stores all information about each particle, i.e.
position, mass, etc...
These are the wall classes.
Referenced by ChuteWithHopper::add_hopper(), compute_walls(), constructor(), Chute::create_bottom(), load_par_ini_file(), print(), read(), read_v1(), read_v2(), Chute::setup_particles_initial_conditions(), ChuteBottom::setup_particles_initial_conditions(), write(), and write_v1().
vector<CWallPeriodic> MD::WallsPeriodic [protected] |
Referenced by Check_and_Duplicate_Periodic_Particle(), do_integration_after_force_computation(), do_integration_before_force_computation(), load_par_ini_file(), print(), read(), read_v1(), read_v2(), Chute::setup_particles_initial_conditions(), ChuteBottom::setup_particles_initial_conditions(), write(), and write_v1().
string MD::xballs_additional_arguments [private] |
Referenced by constructor().
int MD::xballs_cmode [private] |
Referenced by constructor().
Mdouble MD::xballs_scale [private] |
Referenced by constructor().
Mdouble MD::xballs_vscale [private] |
Referenced by constructor().
Referenced by constructor(), output_xballs_data(), print(), read(), read_dim_from_data_file(), read_next_from_data_file(), read_v1(), read_v2(), write(), and write_v1().
These store the size of the domain, assume walls at the ends.
Referenced by constructor(), output_xballs_data(), print(), read(), read_dim_from_data_file(), read_next_from_data_file(), read_v1(), read_v2(), write(), and write_v1().
Referenced by constructor(), output_xballs_data(), print(), read(), read_dim_from_data_file(), read_next_from_data_file(), read_v1(), read_v2(), write(), and write_v1().
Referenced by constructor(), output_xballs_data(), print(), read(), read_dim_from_data_file(), read_next_from_data_file(), read_v1(), read_v2(), write(), and write_v1().
Referenced by constructor(), output_xballs_data(), print(), read(), read_dim_from_data_file(), read_next_from_data_file(), read_v1(), read_v2(), write(), and write_v1().
Referenced by constructor(), output_xballs_data(), print(), read(), read_dim_from_data_file(), read_next_from_data_file(), read_v1(), read_v2(), write(), and write_v1().
1.7.6.1