ChuteBottom.h

Go to the documentation of this file.

//Copyright (c) 2013-2014, The MercuryDPM Developers Team. All rights reserved.
//For the list of developers, see <http://www.MercuryDPM.org/Team>.
//
//Redistribution and use in source and binary forms, with or without
//modification, are permitted provided that the following conditions are met:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name MercuryDPM nor the
//    names of its contributors may be used to endorse or promote products
//    derived from this software without specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
//ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
//WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
//DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
//LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
//ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
//(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
//SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#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 (std::vector<BaseParticle*>::iterator it=getParticleHandler().begin(); it!=getParticleHandler().end(); it++) {
            height = std::max(height,(*it)->get_Position().Z);
        }
        
        std::cout << "Thickness" << thickness << std::endl;
        //now cut a slice of width 2*MaxInflowParticleRadius
        for (std::vector<BaseParticle*>::iterator it=getParticleHandler().begin(); it!=getParticleHandler().end();)
        {
            if ((*it)->get_Position().Z < height - (1.0+thickness)*MaxInflowParticleRadius || (*it)->get_Position().Z > height - MaxInflowParticleRadius) {
                //delete particles outside the given range
                //*it = getParticleHandler().back();
                //getParticleHandler().removeLastObject(); 
                getParticleHandler().removeObject((*it)->get_Index());
            } else {
                //fix the remaining particles on the floor

                //(*it)->get_Position().Z -= height - MaxInflowParticleRadius;
                (*it)->move(Vec3D(0.0,0.0,MaxInflowParticleRadius-height));
                
                (*it)->fixParticle();
                it++;
            }
        }
        
        //copy the rough bottom over
        ChuteParticleHandler.setStorageCapacity(getParticleHandler().getNumberOfObjects());
        std::cout << "Chute bottom finished, consisting of " << getParticleHandler().getNumberOfObjects() << " Particles" << std::endl; 
        
        ChuteParticleHandler = getParticleHandler();
    }

    void setup_particles_initial_conditions() { 
        

        getParticleHandler().setStorageCapacity((int)std::min(3.0*get_xmax()*get_ymax()*get_zmax()/cubic(2.0*get_InflowParticleRadius()),1e6));


        create_bottom();
        
        getWallHandler().clear();
        getBoundaryHandler().clear();
        if (periodicbottom) {
            InfiniteWall w0;
            w0.set(Vec3D(0.0, 0.0, -1.0), -get_zmin()+get_InflowParticleRadius());
            getWallHandler().copyAndAddWall(w0);
            PeriodicBoundary b0;
            b0.set(Vec3D( 1.0, 0.0, 0.0), get_xmin(), get_xmax());
            getBoundaryHandler().copyAndAddObject(b0);
            b0.set(Vec3D( 0.0, 1.0, 0.0), get_ymin(), get_ymax());
            getBoundaryHandler().copyAndAddObject(b0);
        } else {
            InfiniteWall w0;
            w0.set(Vec3D(0.0, 0.0, -1.0), -get_zmin()+get_InflowParticleRadius());
            getWallHandler().copyAndAddWall(w0);
            w0.set(Vec3D(-1.0, 0.0, 0.0), -get_xmin());
            getWallHandler().copyAndAddWall(w0);
            w0.set(Vec3D( 1.0, 0.0, 0.0),  get_xmax());
            getWallHandler().copyAndAddWall(w0);
            w0.set(Vec3D( 0.0,-1.0, 0.0), -get_ymin());
            getWallHandler().copyAndAddWall(w0);
            w0.set(Vec3D( 0.0, 1.0, 0.0),  get_ymax());
            getWallHandler().copyAndAddWall(w0);
        }
        
        //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);

            //The position components are first stored in a Vec3D, because if you pass them directly into set_Position the compiler is allowed to change the order in which the numbers are generated
            Vec3D position;
            position.X=random.get_RN(  P0.get_Radius(), get_xmax()-P0.get_Radius());
            position.Y=random.get_RN(  P0.get_Radius(), get_ymax()-P0.get_Radius());
            position.Z=random.get_RN(2*P0.get_Radius(), get_zmax()-P0.get_Radius());
            P0.set_Position(position);
            P0.set_Velocity(Vec3D(0.0,0.0,0.0));

            if (IsInsertable(P0)) 
            {
                failed = 0; 
                /*if (Particles.capacity()==Particles.size()) 
                {
                    
                    std::cerr << "Error in creating bottom; exceeded capacity: " << Particles.size() << std::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(getParticleHandler().getNumberOfObjects());
        std::cout << "N=" << getParticleHandler().getNumberOfObjects() << std::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.*std::min(get_FixedParticleRadius(),get_MinInflowParticleRadius());
        Mdouble maxCell = 2.*std::max(get_FixedParticleRadius(),get_MaxInflowParticleRadius());
        if ((minCell==maxCell)|(minCell==0.)) setHGridMaxLevels(1);
        else setHGridMaxLevels(2);
//      set_HGRID_cell_to_cell_ratio (1.0000000001*maxCell/minCell);
        //optimize number of buckets
        set_HGRID_num_buckets_to_power(getParticleHandler().getNumberOfObjects()*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 {std::cerr<<"Error: thickness " << new_ << " negative"<<std::endl; exit(-1);}
    }
    Mdouble get_periodicbottom(){return periodicbottom;}
    void set_periodicbottom(bool new_){periodicbottom=new_;}
    
private:
    Mdouble thickness; 
    bool periodicbottom; 
};
#endif