MD.h

Go to the documentation of this file.

#ifndef MD_H
#define MD_H
#define SIGN(x) (x>0?(1):(-1))
#define MAX(x,y) (x>y?(x):(y))
#define MIN(x,y) (x<y?(x):(y))

#define UNUSED  __attribute__ ((__unused__))

//This is a header with some extra standard maths function and constants that are required by the code
#include "ExtendedMath.h"

//The vector class contains a 3D vector class.
#include "Vector.h"
#include <string.h>
//This class defines the particle handler
#include "ParticleHandler.h"
// This class defines walls and periodic walls
#include "CWall.h"
//This is the class that defines the std_save routines
#include "STD_save.h"
#include <time.h>

//This is only used to change the file permission of the xball script create, at some point this code may be moved from this file to a different file.
#include <sys/types.h>
#include <sys/stat.h>

#include "RNG.h"


using namespace std;




class MD : public STD_save {
public:
        void constructor();
        
        MD()
        {
                constructor();
                #ifdef CONSTUCTOR_OUTPUT
                        cerr << "MD() finished"<<endl;
                #endif
        }
        
        MD(STD_save& other) : STD_save(other) {
                constructor();
                #ifdef CONSTUCTOR_OUTPUT
                        cerr << "MD(STD_save& other) finished " << endl;
                #endif          
        }
        
        virtual ~MD() {};
        
        void info(){print(cout);}
        
        void solve();

        void solve(unsigned int argc, char *argv[]) {
                readArguments(argc, argv);
                solve();
        };

        void solveWithMDCLR();
        
        Mdouble get_t(){return t;}
        void set_t(Mdouble new_t){t=new_t;}
        
        void set_NWall(int N){if (N>=0) {Walls.resize(N);} else {cerr << "Error in set_NWall" << endl; exit(-1);}}
        int get_NWall(){return Walls.size();}
        
        //~ ///Allows the number of Species to be changed
        //~ void set_NSpecies(int N){
                //~ if (N>=0) {
                        //~ Species.resize(N); 
                        //~ for (int i=0; i<N; i++) {
                                //~ Species[i].MixedSpecies.resize(i); 
                                //~ cout << "S" << Species[i].MixedSpecies.size() << endl;
                        //~ }
                //~ } else {cerr << "Error in set_NSpecies" << endl; exit(-1);}
        //~ }
        int get_NSpecies(){return Species.size();}
        
        vector<CWall>& get_Walls(void){return Walls;}
        CWall& get_Walls(int i){return Walls[i];}
        
        vector<CWallPeriodic>& get_WallsPeriodic(void){return WallsPeriodic;}
        CWallPeriodic& get_WallsPeriodic(int i){return WallsPeriodic[i];}
        
        vector<CSpecies>& get_Species(void){return Species;}
        CSpecies* get_Species(int i) {return &Species[i];}
        
        CSpecies* get_MixedSpecies(int i, int j) {
                if (i>j) return &Species[i].MixedSpecies[j];
                else return &Species[j].MixedSpecies[i];
        };
        void set_MixedSpecies(int i, int j, CSpecies& S) {
                if (i>j) Species[i].MixedSpecies[j] = S;
                else Species[j].MixedSpecies[i] = S;
        };
        
        void set_NWallPeriodic(int N){if (N>=0) {WallsPeriodic.resize(N);} else {cerr << "Error in set_NWallPeriodic" << endl; exit(-1);}}
        int get_NWallPeriodic(){return WallsPeriodic.size();}
        
        void set_tmax(Mdouble new_tmax){if (new_tmax>=0){tmax=new_tmax;} else { cerr << "Error in set_tmax" << endl; exit(-1); }}
        Mdouble get_tmax(){return tmax;}
        

        ParticleHandler& get_ParticleHandler() { return  particleHandler;}

        
        void set_savecount       (int new_){if (new_>0) {set_save_count_all(new_);}  else {cerr << "Error in set_savecount (set_savecount("<<new_<<"))"<< endl; exit(-1);}}
        void set_save_count_all  (int new_){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 set_save_count_data (int new_){if (new_>0) {save_count_data =new_;} else {cerr << "Error in set_save_count_data, (set_save_count_data ("<<new_<<"))"<< endl; exit(-1);}}
        void set_save_count_ene  (int new_){if (new_>0) {save_count_ene  =new_;} else {cerr << "Error in set_save_count_ene, (set_save_count_ene  ("<<new_<<"))"<< endl; exit(-1);}}
        void set_save_count_stat (int new_){if (new_>0) {save_count_stat =new_;} else {cerr << "Error in set_save_count_stat, (set_save_count_stat ("<<new_<<"))"<< endl; exit(-1);}}
        void set_save_count_fstat(int new_){if (new_>0) {save_count_fstat=new_;} else {cerr << "Error in set_save_count_fstat, (set_save_count_fstat("<<new_<<"))"<< endl; exit(-1);}}
        int get_savecount  (){return get_save_count_data ();}
        int get_save_count (){return get_save_count_data ();}
        int get_save_count_data (){return save_count_data;}
        int get_save_count_ene  (){return save_count_ene;}
        int get_save_count_stat (){return save_count_stat;}
        int get_save_count_fstat(){return save_count_fstat;}
        
        void set_do_stat_always (bool new_) {do_stat_always=new_;}
        void set_number_of_saves      (Mdouble N){set_number_of_saves_all(N);}
        void set_number_of_saves_all  (Mdouble N){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 set_number_of_saves_data (Mdouble N){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 set_number_of_saves_ene  (Mdouble N){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 set_number_of_saves_stat (Mdouble N){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 set_number_of_saves_fstat(Mdouble N){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 set_plastic_k1_k2max_kc_depth(Mdouble k1_, Mdouble k2max_, Mdouble kc_, Mdouble depth_, unsigned int indSpecies = 0)
        {if (indSpecies<Species.size()) {Species[indSpecies].set_plastic_k1_k2max_kc_depth(k1_, k2max_, kc_, depth_);} else {cerr << "Error in set_k: species does not exist"; exit(-1);} }             
        void set_k1(Mdouble new_, unsigned int indSpecies = 0) 
        {if (indSpecies<Species.size()) {Species[indSpecies].set_k1(new_);} else {cerr << "Error in set_k: species does not exist"; exit(-1);} }
        void set_k2max(Mdouble new_, unsigned int indSpecies = 0) 
        {if (indSpecies<Species.size()) {Species[indSpecies].set_k2max(new_);} else {cerr << "Error in set_k: species does not exist"; exit(-1);} }
        void set_kc(Mdouble new_, unsigned int indSpecies = 0) 
        {if (indSpecies<Species.size()) {Species[indSpecies].set_kc(new_);} else {cerr << "Error in set_k: species does not exist"; exit(-1);} }
        void set_depth(Mdouble new_, unsigned int indSpecies = 0) 
        {if (indSpecies<Species.size()) {Species[indSpecies].set_depth(new_);} else {cerr << "Error in set_k: species does not exist"; exit(-1);} }
        Mdouble get_k1(unsigned int indSpecies = 0) 
        {return Species[indSpecies].get_k1();}
        Mdouble get_k2max(unsigned int indSpecies = 0) 
        {return Species[indSpecies].get_k2max();}
        Mdouble get_kc(unsigned int indSpecies = 0) 
        {return Species[indSpecies].get_kc();}
        Mdouble get_depth(unsigned int indSpecies = 0) 
        {return Species[indSpecies].get_depth();}
        Mdouble get_plastic_dt(Mdouble mass, unsigned int indSpecies = 0)
        {return Species[indSpecies].get_plastic_dt(mass);}
        
        void set_k(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_k(new_);} else {cerr << "Error in set_k: species does not exist"; exit(-1);} }
        Mdouble get_k(int indSpecies = 0){return Species[indSpecies].get_k();}
        void set_kt(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_kt(new_);} else {cerr << "Error in set_kt: species does not exist"; exit(-1);}}
        Mdouble get_kt(int indSpecies = 0){return Species[indSpecies].get_kt();}
        void set_krolling(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_krolling(new_);} else {cerr << "Error in set_krolling: species does not exist"; exit(-1);}}
        Mdouble get_krolling(int indSpecies = 0){return Species[indSpecies].get_krolling();}
        void set_ktorsion(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_ktorsion(new_);} else {cerr << "Error in set_ktorsion: species does not exist"; exit(-1);}}
        Mdouble get_ktorsion(int indSpecies = 0){return Species[indSpecies].get_ktorsion();}
        void set_rho(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_rho(new_);} else {cerr << "Error in set_rho: species does not exist"; exit(-1);}}
        Mdouble get_rho(int indSpecies = 0){return Species[indSpecies].get_rho();}
        void set_dispt(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_dispt(new_);} else {cerr << "Error in set_dispt: species does not exist"; exit(-1);}}
        Mdouble get_dispt(unsigned int indSpecies = 0){return Species[indSpecies].get_dispt();}
        void set_disprolling(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_disprolling(new_);} else {cerr << "Error in set_disprolling: species does not exist"; exit(-1);}}
        Mdouble get_disprolling(unsigned int indSpecies = 0){return Species[indSpecies].get_disprolling();}
        void set_disptorsion(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_disptorsion(new_);} else {cerr << "Error in set_disptorsion: species does not exist"; exit(-1);}}
        Mdouble get_disptorsion(unsigned int indSpecies = 0){return Species[indSpecies].get_disptorsion();}

        void set_disp(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_dissipation(new_);} else {cerr << "Error in set_dissipation: species does not exist"; exit(-1);}}
        Mdouble get_disp(unsigned int indSpecies = 0){return Species[indSpecies].get_dissipation();}
        void set_dissipation(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_dissipation(new_);} else {cerr << "Error in set_dissipation: species does not exist"; exit(-1);}}
        Mdouble get_dissipation(unsigned int indSpecies = 0){return Species[indSpecies].get_dissipation();}
        void set_mu(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_mu(new_);} else {cerr << "Error in set_mu: species does not exist"; exit(-1);}}
        Mdouble get_mu(unsigned int indSpecies = 0){return Species[indSpecies].get_mu();}
        void set_murolling(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_murolling(new_);} else {cerr << "Error in set_murolling: species does not exist"; exit(-1);}}
        Mdouble get_murolling(unsigned int indSpecies = 0){return Species[indSpecies].get_murolling();}
        void set_mutorsion(Mdouble new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_mutorsion(new_);} else {cerr << "Error in set_mutorsion: species does not exist"; exit(-1);}}
        Mdouble get_mutorsion(unsigned int indSpecies = 0){return Species[indSpecies].get_mutorsion();}

        void set_dim_particle(int new_, unsigned int indSpecies = 0){if (indSpecies<Species.size()) {Species[indSpecies].set_dim_particle(new_);} else {cerr << "Error in set_dim_particle: species does not exist"; exit(-1);}}
        int get_dim_particle(unsigned int indSpecies = 0){return Species[indSpecies].get_dim_particle();}
        int get_save_data_data(){return save_data_data;}
        int get_save_data_ene(){return save_data_ene;}
        int get_save_data_fstat(){return save_data_fstat;}
        int get_save_data_stat(){return save_data_stat;}
        int get_do_stat_always(){return do_stat_always;}
        
        void set_k_and_restitution_coefficient(Mdouble k_, Mdouble eps, Mdouble mass, unsigned int indSpecies = 0){
                if (indSpecies<Species.size()) {Species[indSpecies].set_k_and_restitution_coefficient(k_, eps, mass);} else {cerr << "Error in set_k_and_restitution_coefficient: species does not exist"; exit(-1);}
        }
        void set_collision_time_and_restitution_coefficient(Mdouble tc, Mdouble eps, Mdouble mass, unsigned int indSpecies = 0){
                if (indSpecies<Species.size()) {Species[indSpecies].set_collision_time_and_restitution_coefficient(tc, eps, mass);} else {cerr << "Error in set_collision_time_and_restitution_coefficient: species does not exist"; exit(-1);}
        }
        
        void set_collision_time_and_restitution_coefficient(Mdouble tc, Mdouble eps, Mdouble mass1, Mdouble mass2,unsigned int indSpecies =0)
                {
                                if (indSpecies<Species.size()) {Species[indSpecies].set_collision_time_and_restitution_coefficient(tc, eps, mass1, mass2);} else {cerr << "Error in set_collision_time_and_restitution_coefficient: species does not exist"; exit(-1);}
                }

        void set_collision_time_and_normal_and_tangential_restitution_coefficient(Mdouble tc, Mdouble eps, Mdouble beta, Mdouble mass1, Mdouble mass2,unsigned int indSpecies =0) {
                if (indSpecies<Species.size()) {Species[indSpecies].set_collision_time_and_normal_and_tangential_restitution_coefficient(tc, eps, beta, mass1, mass2);} else {cerr << "Error in set_collision_time_and_restitution_coefficient: species does not exist"; exit(-1);}
        }

        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) {
                if (indSpecies<Species.size()) {Species[indSpecies].set_collision_time_and_normal_and_tangential_restitution_coefficient_nodispt(tc, eps, beta, mass1, mass2);} else {cerr << "Error in set_collision_time_and_restitution_coefficient: species does not exist"; exit(-1);}

        }

        Mdouble get_collision_time(Mdouble mass, unsigned int indSpecies = 0){return Species[indSpecies].get_collision_time(mass);}
        Mdouble get_restitution_coefficient(Mdouble mass, unsigned int indSpecies = 0){return Species[indSpecies].get_restitution_coefficient(mass);}
        
        Mdouble get_xmin(){return xmin;}
        Mdouble get_xmax(){return xmax;}
        Mdouble get_ymin(){return ymin;}        
        Mdouble get_ymax(){return ymax;}
        Mdouble get_zmin(){return zmin;}        
        Mdouble get_zmax(){return zmax;}
        

        void set_xmin(Mdouble new_xmin){if (new_xmin<=xmax){xmin=new_xmin;} else { cerr << "Error in set_xmin(" << new_xmin << "): xmax=" << xmax << endl; exit(-1); }}

        void set_ymin(Mdouble new_ymin){if (new_ymin<=ymax){ymin=new_ymin;} else { cerr << "Error in set_ymin(" << new_ymin << "): ymax=" << ymax << endl; exit(-1); }}


        void set_zmin(Mdouble new_zmin){if (new_zmin<=zmax){zmin=new_zmin;} else { cerr << "Error in set_zmin(" << new_zmin << "): zmax=" << zmax << endl; exit(-1); }}


        void set_xmax(Mdouble new_xmax){if (new_xmax>=xmin){xmax=new_xmax;} else { cerr << "Error in set_xmax(" << new_xmax << "): xmin=" << xmin << endl; exit(-1); }}

        void set_ymax(Mdouble new_ymax){if (new_ymax>ymin){ymax=new_ymax;} else { cerr << "Error in set_ymax(" << new_ymax << "): ymin=" << ymin << endl; exit(-1); }}

        void set_zmax(Mdouble new_zmax){if (new_zmax>ymin){zmax=new_zmax;} else { cerr << "Error in set_zmax(" << new_zmax << "): zmin=" << zmin << endl; exit(-1); }}
        
        void set_dt(Mdouble new_dt){if (dt>=0.0){dt=new_dt;} else { cerr << "Error in set_dt" << endl; exit(-1); }}
        Mdouble get_dt(){return dt;}
        
        void set_name(const char* name){problem_name.str(""); problem_name << name;}
        
        void set_xballs_colour_mode(int new_cmode){xballs_cmode=new_cmode;}
        void set_xballs_cmode(int new_cmode){xballs_cmode=new_cmode;}
        int get_xballs_cmode(){return xballs_cmode;}
        
        void set_xballs_vector_scale(double new_vscale){xballs_vscale=new_vscale;}
        double get_xballs_vscale(){return xballs_vscale;}
        
        void set_xballs_additional_arguments(string new_){xballs_additional_arguments=new_.c_str();}
        string get_xballs_additional_arguments(){return xballs_additional_arguments;}
        
        void set_xballs_scale(Mdouble new_scale){xballs_scale=new_scale;}
        double get_xballs_scale(){return xballs_scale;}
        
        void set_gravity(Vec3D new_gravity){gravity = new_gravity;}
        Vec3D get_gravity(){return gravity;}
        
        void set_dim(int new_dim){ if (new_dim>=1 && new_dim<=3) dim = new_dim; else { cerr << "Error in set_dim" << endl; exit(-1); }}
        int get_dim(){return dim;}
        
        int get_restart_version(){return restart_version;}
        void set_restart_version(int new_){restart_version=new_;}
        
        bool get_restarted(){return restarted;}

        Mdouble get_max_radius(){return max_radius;}
        
        void set_restarted(bool new_){
                restarted=new_;
                set_append(new_);
        }

        bool get_append(){return append;}
        void set_append(bool new_){append=new_;}


        Mdouble get_ene_ela(){return ene_ela;}
        void set_ene_ela(Mdouble new_){ene_ela=new_;}
        void add_ene_ela(Mdouble new_){ene_ela+=new_;}

        //access function to MD::Hertzian
        bool get_Hertzian(){return Hertzian;}
        void set_Hertzian(bool new_){Hertzian=new_;
                if (Hertzian) cout << "Hertzian law activated" << endl;
                else cout << "Hertzian law deactivated" << endl;
        }

        void Remove_Particle(int IP);
        
        // Some statistical output that might be of interest when setting particles' parameters
        
        Mdouble get_Mass_from_Radius(Mdouble radius, int indSpecies=0) {
                Particle P;
                P.set_Radius(radius);
                P.set_IndSpecies(indSpecies);
                P.compute_particle_mass(Species);
                return P.get_Mass();
        }
        
        Mdouble get_maximum_velocity(Particle &P){return Species[P.get_IndSpecies()].get_maximum_velocity(P.get_Radius(),P.get_Mass());}
        
        virtual Particle* get_SmallestParticle() {return particleHandler.get_SmallestParticle();}
        virtual Particle* get_LargestParticle() {return particleHandler.get_LargestParticle();}
        virtual void removeParticle(int iP)
        {
                HGRID_RemoveParticleFromHgrid(particleHandler.get_Particle(iP));
                particleHandler.removeParticle(iP);
        }
        
        Mdouble get_maximum_velocity(){
                Particle *P = get_SmallestParticle();
                return P->get_Radius() * sqrt(Species[P->get_IndSpecies()].k/(.5*P->get_Mass()));
        }
        
        void set_dt_by_mass(Mdouble mass){dt = .02 * get_collision_time(mass);}
        
        void set_dt(Particle &P){
                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 set_dt(){
                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]);
        }
        
        virtual void setup_particles_initial_conditions();
        
        virtual void create_xballs_script();

        virtual double getInfo(Particle& P UNUSED) {return 0;}
        
        virtual void save_restart_data ();
        
        int load_restart_data ();
        int load_restart_data (string filename);
        
        void statistics_from_restart_data(const char* name);
        
        virtual void write(std::ostream& os);
        virtual void read(std::istream& is);
        virtual void write_v1(std::ostream& os);
        virtual void read_v1(std::istream& is);
        virtual void read_v2(std::istream& is);
        
        bool load_from_data_file (const char* filename, unsigned int format=0);
        bool load_par_ini_file (const char* filename);
        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);
        
        friend inline std::ostream& operator<<(std::ostream& os, MD &md) {
                md.print(os);
                return os;
        }
        
        void add_Species(CSpecies& S);
        
        void add_Species(void) { add_Species(Species[0]); };
                
        void set_format(int new_) {format = new_;};
        
        int readArguments(unsigned int argc, char *argv[]);

        virtual int readNextArgument(unsigned int& i, unsigned int argc, char *argv[]);
        
        //This is a random generator, often used by the initial conditions etc...
        RNG random;

//functions that should only be used in the class definitions 
protected:
        
        virtual void compute_all_forces();

        virtual void compute_internal_forces(Particle* i);

        virtual void compute_internal_forces(Particle* P1, Particle* P2);
        
        void compute_plastic_internal_forces(Particle* P1, Particle* P2);
        
        virtual void compute_external_forces(Particle* PI);
        
        virtual void compute_walls(Particle* PI);
        
        virtual void actions_before_time_loop(){
                //automatically sets dt if dt is not specified by the user
                if (!dt) set_dt(); 
        };
        
        virtual void HGRID_actions_before_time_loop(){};
        
        virtual void HGRID_actions_before_time_step(){};
        
        virtual void HGRID_InsertParticleToHgrid(Particle *obj UNUSED){};
        virtual void HGRID_UpdateParticleInHgrid(Particle *obj UNUSED){};
        virtual void HGRID_RemoveParticleFromHgrid(Particle *obj UNUSED){};
        virtual bool get_HGRID_UpdateEachTimeStep(){return true;};
        
        virtual void actions_before_time_step(){};
        
        virtual void actions_after_solve(){};

        virtual void actions_after_time_step(){};       
        
        virtual void output_xballs_data();
        virtual void output_xballs_data_particle(int i);
        
        virtual void start_ene();
        virtual void fstat_header();
        virtual void output_ene();
        
        virtual void initialize_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);
        
        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);
        
        virtual void InitBroadPhase(){};
        
        virtual void broad_phase(Particle* i)
        {
                for (vector<Particle*>::iterator it = particleHandler.begin(); (*it)!=i; ++it)
                {
                        compute_internal_forces(i,*it);
                }
        }
        
        void set_FixedParticles(unsigned int n){for (unsigned int i=0; i<std::min((unsigned int)particleHandler.get_NumberOfParticles(),n); i++) particleHandler.get_Particle(i)->fixParticle();}
        
        void initialize_tangential_springs();
                
        void compute_particle_masses() {for (unsigned int i=0;i<particleHandler.get_NumberOfParticles();i++) particleHandler.get_Particle(i)->compute_particle_mass(Species);}
        
        virtual void cout_time() {
                cout << "\rt=" << setprecision(3) <<fixed<< left << setw(6) << t 
                        << ", tmax=" << setprecision(3) << left << setw(6) << tmax << " \r";
                cout.flush();
        }
        
        virtual bool continue_solve() {return true;}
        
        void reset_DeltaMax(){for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) (*it)->get_DeltaMaxs().resize(0);}
        void reset_TangentialSprings(){for (vector<Particle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) (*it)->get_TangentialSprings().resize(0);}


        //virtual int Check_and_Duplicate_Periodic_Particle(int i, int nWallPeriodic);

        //vector<Particle> Particles;
        
        vector<CWall> Walls;
        vector<CWallPeriodic> WallsPeriodic;

        vector<CSpecies> Species;

private:

        int dim;

        Vec3D gravity; 
        
        Mdouble max_radius;
        
        Mdouble xmin,xmax,ymin,ymax,zmin,zmax;
        
        Mdouble dt;
        Mdouble tmax;
        
        int save_count_data; //<determines how many time steps are skipped before the next data file is written
        int save_count_ene;  //<determines how many time steps are skipped before the next ene file is written
        int save_count_stat; //<determines how many time steps are skipped before the next fstat file is written
        int save_count_fstat;//<determines how many time steps are skipped before the next stat file is written
        bool save_data_data;//<determines if the current timestep is written into the data file
        bool save_data_ene;//<determines if the current timestep is written into the ene file
        bool save_data_fstat;//<determines if the current timestep is written into the fstat file
        bool save_data_stat;//<determines if the current timestep is written into the stat file
        bool do_stat_always;
        

        Mdouble t;

        Mdouble ene_ela; //<stores the potential energy in the elastic springs (global, because it has to be calculated in the force loop

        
        //This is the private data that is only used by the xballs output

        int xballs_cmode; //< sets the xballs argument cmode (see xballs.txt)
        Mdouble xballs_vscale; //< sets the xballs argument vscale (see xballs.txt)
        Mdouble xballs_scale; //< sets the xballs argument scale (see xballs.txt)
        string xballs_additional_arguments; //< string where additional xballs argument can be specified (see xballs.txt)
        int format; //< sets the data file argument format (i.e. how many columns the data file has)
        int restart_version; //<to read new and old restart data
        bool restarted; //<to read new and old restart data
        bool Hertzian; //<make force law Hertzian
        int data_FixedParticles; //<determines how many particles are fixed when reading in data files

        
private:

        ParticleHandler particleHandler;
        
        void Remove_Duplicate_Periodic_Particles();
        virtual int Check_and_Duplicate_Periodic_Particle(Particle* i, int nWallPeriodic);
        int Check_and_Duplicate_Periodic_Particles();


        bool append;
        bool rotation;
        int PeriodicCreated;
        int save_restart_data_counter;

        


};

#endif