HGRID_base.h

Go to the documentation of this file.

#ifndef HGRID_H
#define HGRID_H

#include <ostream>
#include <iostream>
#include <stdlib.h>
using std::cout;
using std::cerr;

#include "CParticle.h"
#include "MD.h"



class Cell 
{
public:
Cell () { }
Cell(int px, int py, int pz, int pl) { x = px; y = py; z = pz; l = pl; }

int x, y, z, l;


friend inline std::ostream& operator<<(std::ostream& os, const Cell &cell)
        {
        os << cell.x << ' ' << cell.y << ' ' << cell.z << ' ' << cell.l;
        return os;
        }
};


class HGrid
{
public:
        
        HGrid() {
                objectBucket = NULL; 
                bucketIsChecked = NULL;
                //minCellPos = NULL; 
                //maxCellPos = NULL; 
                inv_size = NULL; 
        }
        
        /*void print(std::ostream& os) {
                os << "HGRID" << endl;
                for (int i=0; i<NUM_BUCKETS; i++)
                {
                        if (objectBucket[i]!=-1) {
                                os << i << " ";
                                for (int p=objectBucket[i]; p!=-1; p=Particles[p]->pNextObject)
                                        os << p << " ";
                                os << endl;
                        }
                }
        }*/
        
        
        //----- parameters ------
        int NUM_BUCKETS; 
        
        int HGRID_MAX_LEVELS; 
        
        Mdouble MIN_CELL_SIZE; 
        
        Mdouble SPHERE_TO_CELL_RATIO; 
                
        Mdouble CELL_TO_CELL_RATIO; 
        
        //----- internal variables ----------
        int occupiedLevelsMask; 
        
        Particle **objectBucket; 
        
        bool *bucketIsChecked; 
        
        //Cell *minCellPos; 
        
        //Cell *maxCellPos; 
        
        Mdouble *inv_size; 
        
        //-------------------------------------------------
        HGrid(int num_buckets, int hgrid_max_levels, Mdouble min_cell_size, Mdouble sphere_to_cell_ratio, Mdouble cell_to_cell_ratio)
        {
                NUM_BUCKETS = num_buckets;
                
                objectBucket = new Particle*[NUM_BUCKETS];
                bucketIsChecked = new bool[NUM_BUCKETS];
                
                HGRID_MAX_LEVELS = hgrid_max_levels;
                
                //minCellPos = new Cell[HGRID_MAX_LEVELS];
                //maxCellPos = new Cell[HGRID_MAX_LEVELS];
                inv_size = new Mdouble[HGRID_MAX_LEVELS];
                
                
                MIN_CELL_SIZE = min_cell_size;
                
                SPHERE_TO_CELL_RATIO = sphere_to_cell_ratio;
                
                CELL_TO_CELL_RATIO = cell_to_cell_ratio;
        }

        ~HGrid() {
                if (objectBucket) { delete [] objectBucket; objectBucket=NULL; }
                if (bucketIsChecked) { delete [] bucketIsChecked; bucketIsChecked=NULL; }
                //if (minCellPos) { delete [] minCellPos; minCellPos=NULL; }
                //if (maxCellPos) { delete [] maxCellPos; maxCellPos=NULL; }
                if (inv_size) { delete [] inv_size; inv_size=NULL; } 
        }
        
        //-------------------------------------------------
        void Initialize_inv_size()
        {
                Mdouble size = MIN_CELL_SIZE;
                
                for (int level = 0; level<HGRID_MAX_LEVELS; level++)
                {
                        // initialize inv_size
                        this->inv_size[level] = 1.0f/size;
                                                
                        size *= CELL_TO_CELL_RATIO;
                }
        }
        
        //-------------------------------------------------
        /*void InitializeMaxMinCellPos(Mdouble wx0, Mdouble wy0, Mdouble wz0, Mdouble wx1, Mdouble wy1, Mdouble wz1)
        {
                Mdouble inv_size;
                                        
                for (int level = 0; level<HGRID_MAX_LEVELS; level++)
                {
                        inv_size = this->inv_size[level];
                                        
                        this->minCellPos[level] = Cell((int)(wx0 * inv_size), (int)(wy0 * inv_size), (int)(wz0 * inv_size), level);
                        this->maxCellPos[level] = Cell((int)(wx1 * inv_size), (int)(wy1 * inv_size), (int)(wz1 * inv_size), level);
                }
        
        }*/
        
        void InsertParticleToHgrid(Particle *obj)
        {
                // Find lowest level where object fully fits inside cell, taking RATIO into account
                int level;
                Mdouble size = MIN_CELL_SIZE;
                Mdouble diameter = obj->get_Radius() * 2.0f;
                
                for (level = 0; size  < SPHERE_TO_CELL_RATIO*diameter; level++)
                {
                        size *= CELL_TO_CELL_RATIO;
                }

                if (level >= HGRID_MAX_LEVELS) 
                { 
                        cerr << "ATTENTION !!! Object is larger (d=" << diameter << ") than largest grid cell (" << size << ")!" << std::endl; 
                        cerr << "Please lower minimum cell size (" << MIN_CELL_SIZE << ") for HGRID and rerun" << std::endl;
                        //cerr << "HGRID_MAX_LEVELS" << HGRID_MAX_LEVELS << "CELL_TO_CELL_RATIO" << CELL_TO_CELL_RATIO << endl;
                        exit(-1);
                }
                        
                obj->set_HGRID_Level(level);
                        
                // indicate level is in use - not levels with no particles no collision detection is performed
                this->occupiedLevelsMask |= (1 << level);
        }
        
                                                //-------------------------------------------------
                
        int ComputeHashBucketIndex(Cell cellPos)
        {
                const int h1 = 0x8da6b343; // Large multiplicative constants;
                const int h2 = 0xd8163841; // here arbitrarily chosen primes
                const int h3 = 0xcb1ab31f;
                const int h4 = 0x165667b1;
                        
                int n = h1 * cellPos.x + h2 * cellPos.y + h3 * cellPos.z + h4 * cellPos.l;
                n = n % NUM_BUCKETS;
                        
                if (n < 0) n += NUM_BUCKETS;
                                                                                                                        
                return n;
        }

        void reset_num_buckets(int new_num_buckets){
                NUM_BUCKETS=new_num_buckets;
                objectBucket = (Particle**)realloc(objectBucket, NUM_BUCKETS*sizeof(Particle*));
                bucketIsChecked = (bool*)realloc(bucketIsChecked, NUM_BUCKETS*sizeof(bool));
        }
};
//(-)********** class HGrid *****************************************
        
        
class HGRID_base : public virtual MD
{

public:
        
        HGRID_base()
        {
                constructor();
                #ifdef CONSTUCTOR_OUTPUT
                        cerr << "HGRID_base() finished"<<endl;
                #endif
        }

        ~HGRID_base() { 
                if (grid) { delete grid; grid=NULL; } 
        }
        
        HGRID_base(MD& other) : MD(other)
        {
                constructor();
                #ifdef CONSTUCTOR_OUTPUT
                        cerr << "HGRID_base(MD& other) finished " << endl;
                #endif
        }
        
        void constructor();
        
        void HGRID_actions_before_time_loop(){InitBroadPhase();}
        
        void HGRID_actions_before_time_step();
        
        void set_HGRID_num_buckets(unsigned int new_num_buckets){if (new_num_buckets>0) NUM_BUCKETS=new_num_buckets; else {cerr << "Error in set_HGRID_num_buckets" << endl; exit(-1);}}
        
        void set_HGRID_num_buckets_to_power(){ set_HGRID_num_buckets_to_power(get_ParticleHandler().get_NumberOfParticles());   }
        
        void set_HGRID_num_buckets_to_power(unsigned int N){
                unsigned int NUM_BUCKETS = 32;
                while (NUM_BUCKETS<2*N) NUM_BUCKETS *= 2;
                set_HGRID_num_buckets(NUM_BUCKETS);
        }

        void set_HGRID_max_levels(int new_max_levels){if (new_max_levels>0){HGRID_MAX_LEVELS=new_max_levels;}}
        
        void set_HGRID_sphere_to_cell_ratio(Mdouble new_ratio){if (new_ratio>0){SPHERE_TO_CELL_RATIO=new_ratio;}}
        
        void set_HGRID_cell_to_cell_ratio(Mdouble new_ratio){if (new_ratio>0){CELL_TO_CELL_RATIO=new_ratio;}}
        
        void read(std::istream& is)
        { 
                MD::read(is); 


                //read out the full line first, so if there is an error it does not affect the read of the next line
                string line_string;
                getline(is,line_string);
                stringstream line (stringstream::in | stringstream::out);
                line << line_string;

                if (get_restart_version()==1) {
                line >> NUM_BUCKETS >> HGRID_MAX_LEVELS >> MIN_CELL_SIZE >> SPHERE_TO_CELL_RATIO >> CELL_TO_CELL_RATIO;
                } else {
      string dummy;
        line>> dummy >> NUM_BUCKETS 
        >> dummy >> HGRID_MAX_LEVELS 
        >> dummy >> MIN_CELL_SIZE 
        >> dummy >> SPHERE_TO_CELL_RATIO 
        >> dummy >> CELL_TO_CELL_RATIO;
    }
        }
        
        void write(std::ostream& os)
        { 
                MD::write(os); 
                os      << "NUM_BUCKETS " << NUM_BUCKETS << " "
                        << "HGRID_MAX_LEVELS " << HGRID_MAX_LEVELS << " "
                        << "MIN_CELL_SIZE " << MIN_CELL_SIZE << " "
                        << "SPHERE_TO_CELL_RATIO " << SPHERE_TO_CELL_RATIO  << " "
                        << "CELL_TO_CELL_RATIO " << CELL_TO_CELL_RATIO << endl;
        }
        
  void print(std::ostream& os, bool print_all=false)
  {
                MD::print(os,print_all);
                os << " NUM_BUCKETS:" << NUM_BUCKETS << ", HGRID_MAX_LEVELS:" << HGRID_MAX_LEVELS << ", MIN_CELL_SIZE:" << MIN_CELL_SIZE << ", " << endl
                   << " SPHERE_TO_CELL_RATIO:" << SPHERE_TO_CELL_RATIO << ", CELL_TO_CELL_RATIO:" << CELL_TO_CELL_RATIO << endl;
  }
  
        Mdouble get_CurrentMaxRelativeDisplacement(){return currentMaxRelativeDisplacement;}
        Mdouble get_TotalCurrentMaxRelativeDisplacement(){return totalCurrentMaxRelativeDisplacement;}
        
        void set_UpdateEachTimeStep(bool _new){updateEachTimeStep=_new;}
        bool get_HGRID_UpdateEachTimeStep(){return updateEachTimeStep;}   

        protected:
        
        void InitBroadPhase();
        
        void HGRID_InsertParticleToHgrid(Particle *obj);
                        
        void broad_phase(Particle *i){CheckObjAgainstGrid(grid, i);}
        
        virtual void CheckObjAgainstGrid(HGrid *grid, Particle *obj)=0;
        
        virtual bool TestObject(Particle *pI, Particle *pJ) {
                //PI==PJ check is required because this function is called for all possible close Particle combination (even itself)
                return pI->get_Index()==pJ->get_Index() || GetDistance2(pI->get_Position(),pJ->get_Position())>=sqr(pI->get_Radius()+pJ->get_Radius());
        }


        void HGRID_update_move(Particle * iP, Mdouble move)
        {
                Mdouble currentRelativeDisplacement=move/(grid->MIN_CELL_SIZE*pow(grid->CELL_TO_CELL_RATIO,iP->get_HGRID_Level()));
                if(currentRelativeDisplacement>currentMaxRelativeDisplacement) currentMaxRelativeDisplacement=currentRelativeDisplacement;
        };
        
        void HGRID_actions_before_integration(){currentMaxRelativeDisplacement=0;}
        void HGRID_actions_after_integration(){totalCurrentMaxRelativeDisplacement+=2.0*currentMaxRelativeDisplacement;}
        
        private:
        
        void set_HGRID_min_cell_size(Mdouble new_min_cell_size){if (new_min_cell_size>0){MIN_CELL_SIZE=new_min_cell_size;}}
        void set_CELL_TO_CELL_RATIO(Mdouble new_CELL_TO_CELL_RATIO){if (new_CELL_TO_CELL_RATIO>0){CELL_TO_CELL_RATIO=new_CELL_TO_CELL_RATIO;}}
        
        Mdouble currentMaxRelativeDisplacement, totalCurrentMaxRelativeDisplacement;
        bool updateEachTimeStep;
        int NUM_BUCKETS;
        int HGRID_MAX_LEVELS;
        Mdouble MIN_CELL_SIZE;
        Mdouble SPHERE_TO_CELL_RATIO;
        Mdouble CELL_TO_CELL_RATIO;
        
        
        
        protected:
        
        
        
        HGrid* grid;
        
        void fix_hgrid (){
                set_HGRID_min_cell_size(2.00001*get_SmallestParticle()->get_Radius());
                set_HGRID_max_levels(2);
                set_CELL_TO_CELL_RATIO(get_LargestParticle()->get_Radius()/get_SmallestParticle()->get_Radius());
                set_HGRID_num_buckets_to_power(get_ParticleHandler().get_NumberOfParticles());
        }
};



#endif