ChuteBottom.h

Go to the documentation of this file.

#ifndef CHUTEBOTTOM_H
#define CHUTEBOTTOM_H
#include "Chute.h"

class ChuteBottom : public Chute {

public:
        ChuteBottom() {constructor();}
        
        ChuteBottom(MD& other) :                        MD(other), Chute(other) {constructor();}
        ChuteBottom(HGRID_base& other) :        MD(other), Chute(other) {constructor();}
        ChuteBottom(HGRID_3D& other) :          MD(other), Chute(other) {constructor();}
        ChuteBottom(Chute& other) :                     MD(other), Chute(other) {constructor();}
        
        void constructor(){     
                set_name("roughbottom");
                set_options_fstat(0); //set to 0 for no data creation
                set_options_data(1);
                set_options_restart(1);
                set_thickness(2.4);
                set_periodicbottom(true);
        }

        void make_rough_bottom(ParticleHandler& ChuteParticleHandler) 
        {
                // set all parameters that should be different from the original chute 
                set_ChuteAngle(0.0);
                set_InflowHeight(25.*get_InflowParticleRadius());
                //~ set_InflowHeight(45.*get_InflowParticleRadius()); note: Changing the Inflow height was an attempt to make the bottom density homogeneous, but it did not have the desired effect
                set_RandomizedBottom(1);
                set_FixedParticleRadius(get_InflowParticleRadius());

                set_collision_time_and_restitution_coefficient(get_collision_time()*10.0, 0.2);
                set_mu(0);
                
                set_dt();
                set_dt(get_dt()*10.0);
                set_tmax(get_dt()*2e3);
                //set_number_of_saves(2);
                set_savecount(100);

                //solve
                solve();

                //create_bottom         
                Mdouble height = 0;
                for (vector<Particle*>::iterator it=get_ParticleHandler().begin(); it!=get_ParticleHandler().end(); it++) {
                        height = max(height,(*it)->get_Position().Z);
                }
                
                cout << "Thickness" << thickness << endl;
                //now cut a slice of width 2*MaxInflowParticleRadius
                for (vector<Particle*>::iterator it=get_ParticleHandler().begin(); it!=get_ParticleHandler().end();) {
                        if ((*it)->get_Position().Z < height - (1.0+thickness)*MaxInflowParticleRadius || (*it)->get_Position().Z > height - MaxInflowParticleRadius) {
                                //delete particles outside the given range
                                //*it = get_ParticleHandler().back();
                                //get_ParticleHandler().removeLastParticle(); 
                                get_ParticleHandler().removeParticle((*it)->get_Index());
                        } else {
                                //fix the remaining particles on the floor
                                (*it)->get_Position().Z -= height - MaxInflowParticleRadius;
                                (*it)->fixParticle();
                                it++;
                        }
                }
                
                //copy the rough bottom over
                ChuteParticleHandler.set_StorageCapacity(get_ParticleHandler().get_NumberOfParticles());
                cout << "Chute bottom finished, consisting of " << get_ParticleHandler().get_NumberOfParticles() << " Particles" << endl; 
                
                ChuteParticleHandler = get_ParticleHandler();
        }

        void setup_particles_initial_conditions() {     
                

                get_ParticleHandler().set_StorageCapacity((int)min(3.0*get_xmax()*get_ymax()*get_zmax()/cubic(2.0*get_InflowParticleRadius()),1e6));


                create_bottom();
                
                if (periodicbottom) {
                        Walls.resize(1);
                        Walls.back().set(Vec3D(0.0, 0.0, -1.0), -get_zmin()+get_InflowParticleRadius());
                        WallsPeriodic.resize(2);
                        WallsPeriodic[0].set(Vec3D( 1.0, 0.0, 0.0), get_xmin(), get_xmax());
                        WallsPeriodic[1].set(Vec3D( 0.0, 1.0, 0.0), get_ymin(), get_ymax());
                } else {
                        WallsPeriodic.resize(0);
                        Walls.resize(5);
                        Walls[0].set(Vec3D(0.0, 0.0, -1.0), -get_zmin()+get_InflowParticleRadius());
                        Walls[Walls.size()-4].set(Vec3D(-1.0, 0.0, 0.0), -get_xmin());
                        Walls[Walls.size()-3].set(Vec3D( 1.0, 0.0, 0.0),  get_xmax());
                        Walls[Walls.size()-2].set(Vec3D( 0.0,-1.0, 0.0), -get_ymin());
                        Walls[Walls.size()-1].set(Vec3D( 0.0, 1.0, 0.0),  get_ymax());
                }
                
                //add particles
                HGRID_actions_before_time_loop();
                int failed = 0, max_failed = 500;
                //try max_failed times to find new insertable particle
                while (failed<=max_failed){
                        P0.set_Radius(FixedParticleRadius);
                        P0.compute_particle_mass(Species);

                        P0.get_Position().X = random.get_RN(P0.get_Radius(), get_xmax()-P0.get_Radius());
                        P0.get_Position().Y = random.get_RN(P0.get_Radius(), get_ymax()-P0.get_Radius());
                        P0.get_Position().Z = random.get_RN(2*P0.get_Radius(), get_zmax()-P0.get_Radius());
                        P0.set_Velocity(Vec3D(0.0,0.0,0.0));

                        if (IsInsertable(P0)) 
                        {
                                failed = 0; 
                                /*if (Particles.capacity()==Particles.size()) 
                                {
                                        
                                        cerr << "Error in creating bottom; exceeded capacity: " << Particles.size() << endl; exit(-1);
                                }*/
                                //Particles.push_back (P0); 
                                //Undate the link list so this particle is back in the grid
                                //grid->objectBucket[Particles[Particles.size()-1].bucket] = Particles.size()-1;
                                //grid->objectBucket[Particles[Particles.size()-1].bucket] = &Particles[Particles.size()-1];
                                num_created++;
                        } 
                        else failed++;
                }
                //set_Nmax(get_ParticleHandler().get_NumberOfParticles());
                cout << "N=" << get_ParticleHandler().get_NumberOfParticles() << endl;

                //fix hgrid (there is still an issue when particles are polydispersed)
                //assume 1-2 levels are optimal (which is the case for mono and bidispersed) and set the cell size to min and max 
                // !this is not optimal for polydispersed
                Mdouble minCell = 2.*min(get_FixedParticleRadius(),get_MinInflowParticleRadius());
                Mdouble maxCell = 2.*max(get_FixedParticleRadius(),get_MaxInflowParticleRadius());
                if ((minCell==maxCell)|(minCell==0.)) set_HGRID_max_levels(1);
                else set_HGRID_max_levels(2);
                set_HGRID_cell_to_cell_ratio (1.0000000001*maxCell/minCell);
                //optimize number of buckets
                set_HGRID_num_buckets_to_power(get_ParticleHandler().get_NumberOfParticles()*1.5);
                //end: fix hgrid
                
                //~ info();
        }
                
        void actions_before_time_step() {       };
        
        Mdouble get_thickness(){return thickness;}
        void set_thickness(Mdouble new_){
                if (new_>0.0) thickness=new_;
                else {cerr<<"Error: thickness " << new_ << " negative"<<endl; exit(-1);}
        }
        Mdouble get_periodicbottom(){return periodicbottom;}
        void set_periodicbottom(bool new_){periodicbottom=new_;}
        
private:
        Mdouble thickness; 
        bool periodicbottom; 
};
#endif