box.hpp

Go to the documentation of this file.

#include "vecrand.hpp"

#ifndef BOX_H
#define BOX_H

#include <vector>
#include <boost/shared_ptr.hpp>

#define sptr boost::shared_ptr

using namespace boost; // required for SWIG for some reason

typedef const unsigned int cuint;

inline bool to_buffer(vector<Vec*> arr, double* buffer, size_t sizet) {
    if(sizet < NDIM * arr.size()){return false;};
    for(uint i=0; i < arr.size(); i++)
    for(uint j=0; j < NDIM; j++)
    {
        buffer[i*NDIM + j] = (double) ((*arr[i])[j]);
    }
    return true;
};

class AtomGroup;

class Box {
    public:

        virtual Vec diff(Vec r1, Vec r2)=0;
        virtual flt V()=0;
        virtual ~Box(){};
};

/***********************************************************************
 * Boxes
 */


#ifdef VEC3D
inline Vec vec_mod(Vec r1, Vec r2){
    return Vec(remainder(r1[0], r2[0]), remainder(r1[1], r2[1]), remainder(r1[2], r2[2]));
};
#endif
#ifdef VEC2D
inline Vec vec_mod(Vec r1, Vec r2){
    return Vec(remainder(r1[0], r2[0]), remainder(r1[1], r2[1]));
};
#endif


class InfiniteBox : public Box {
    public:
        Vec diff(Vec r1, Vec r2){return r1-r2;};
        flt V(){return NAN;};
};


class OriginBox : public Box {
    protected:
        Vec boxsize;
    public:
        OriginBox(Vec size) : boxsize(size){};
        Vec diff(Vec r1, Vec r2){
            return vec_mod((r1-r2), boxsize);
        };
        virtual Vec diff(Vec r1, Vec r2, boost::array<int,NDIM> boxes){
            Vec dr = r1 - r2;
            for(uint i=0; i<NDIM; i++) dr[i] -= boxsize[i] * boxes[i];
            return dr;
        };
        virtual boost::array<int,NDIM> box_round(Vec r1, Vec r2){
            boost::array<int,NDIM> boxes;
            Vec dr = r1 - r2;
            for(uint i=0; i<NDIM; i++) boxes[i] = (int) round(dr[i] / boxsize[i]);
            return boxes;
        };
        #ifdef VEC3D
        OriginBox(flt L) : boxsize(L,L,L){};
        flt V(){return boxsize[0] * boxsize[1] * boxsize[2];};
        flt L(){return (boxsize[0] + boxsize[1] + boxsize[2])/3.0;};
        #endif
        #ifdef VEC2D
        OriginBox(flt L) : boxsize(L,L){};
        flt V(){return boxsize[0] * boxsize[1];};
        flt L(){return (boxsize[0] + boxsize[1])/2.0;};
        #endif
        flt resize(flt factor);
        flt resize(flt factor, AtomGroup &atoms);
        flt resize_to(Vec newsize);
        flt resize_to(Vec newsize, AtomGroup &atoms);
        flt resize_to_V(flt newV);
        flt resize_to_V(flt newV, AtomGroup &atoms);
        
        flt resize_to_L(flt newL);
        flt resize_to_L(flt newL, AtomGroup &atoms);
        Vec rand_loc(){
            Vec v = rand_vec_boxed();
            for(uint i=0; i<NDIM; i++){
                v[i] *= boxsize[i];
            }
            return diff(v, Vec::Zero());
        };
        Vec box_shape(){return boxsize;};
        
        void pure_shear_to(flt epsilon);
        void pure_shear_to(flt epsilon, AtomGroup &atoms);
};


class LeesEdwardsBox : public OriginBox {
    protected:
        flt gamma;
    public:
        LeesEdwardsBox(Vec size, flt gamma=0.0) : OriginBox(size), gamma(gamma){};
        LeesEdwardsBox(flt L, flt gamma=0.0) : OriginBox(L), gamma(gamma){};
        Vec diff(Vec r1, Vec r2);
        virtual Vec diff(Vec r1, Vec r2, boost::array<int,NDIM> boxes);
        virtual boost::array<int,NDIM> box_round(Vec r1, Vec r2);

        flt get_gamma(){return gamma;};

        void shear(flt dgamma);
        void shear(flt dgamma, AtomGroup &atoms);
        void shear_to(flt gamma);
        void shear_to(flt gamma, AtomGroup &atoms);
        Vec non_affine(Vec v){
            v[0] -= gamma * v[1];
            return v;
        }

        Vec affine(Vec v){
            v[0] += gamma * v[1];
            return v;
        }
};


class SCBox : public Box {
    protected:
        flt L, R;
    public:
        SCBox(flt L, flt R);
        Vec diff(Vec r1, Vec r2){return r1-r2;};
        flt V();
        Vec dist(Vec r1);
        Vec edge_dist(Vec r1);
        bool inside(Vec r1, flt buffer=0.0);
        Vec rand_loc(flt min_dist_to_wall=0.0);
        flt length(){return L;};
        flt radius(){return R;};
};


struct Atom {
    Vec x;

    Vec v;

    Vec a;

    Vec f;

    flt m;
};

class AtomRef {
    private:
        Atom *ptr;
    public:
        inline AtomRef() : ptr(NULL){};
        inline AtomRef(Atom *a) : ptr(a){};
        inline Atom& operator *() const {return *ptr;};
        inline Atom *operator->() const{ return ptr;}
        inline bool operator==(const AtomRef &other) const {return other.ptr == ptr;};
        inline bool operator==(const Atom* other) const {return other == ptr;};
        inline bool operator!=(const AtomRef &other) const {return other.ptr != ptr;};
        inline bool operator<(const AtomRef &other) const {return ptr < other.ptr;};
        inline bool operator<=(const AtomRef &other) const {return ptr <= other.ptr;};
        inline bool operator>=(const AtomRef &other) const {return ptr >= other.ptr;};
        inline bool operator>(const AtomRef &other) const {return ptr > other.ptr;};
        inline bool is_null(){return ptr==NULL;};
};


class AtomID : public AtomRef {
    private:
        uint num; // note that these are generally only in reference to
                  // a specific AtomGroup
    public:
        inline AtomID() : AtomRef(), num(UINT_MAX){};
        inline AtomID(Atom *a, uint n) : AtomRef(a), num(n){};
        inline uint n() const {return num;};
};

class IDPair {
    private:
        AtomID id1, id2;
    public:
        IDPair() : id1(), id2(){};
        IDPair(AtomID a, AtomID b) : id1(a), id2(b){};
        inline AtomID first() const {return id1;};
        inline AtomID last() const {return id2;};
};

class AtomGroup;

class AtomIter{
    private:
        uint i;
        AtomGroup &g;
    public:
        AtomIter (AtomGroup& g, uint i): i(i), g(g){};
        bool operator!=(const AtomIter& other) const{return i != other.i;};
        Atom& operator* () const;
        inline const AtomIter& operator++(){++i; return *this;};
};

class AtomVec;


class AtomGroup {
    public:
        // access individual atoms
        virtual AtomVec& vec()=0;
        virtual Atom& operator[](cuint n)=0;
        virtual Atom& operator[](cuint n) const=0;
        virtual Atom& get(cuint n){return ((*this)[n]);};
        virtual AtomID get_id(cuint n)=0;
        virtual uint size() const=0;
        virtual AtomIter begin(){return AtomIter(*this, 0);};
        virtual AtomIter end(){return AtomIter(*this, (uint) size());};

        Vec com() const;
        Vec com_force() const;
        Vec com_velocity() const;

        flt mass() const;

        flt kinetic_energy(const Vec originvelocity=Vec::Zero()) const;
        Vec momentum() const;
        flt gyradius() const;
        #ifdef VEC3D
        Vec torque(const Vec loc) const;
        Vec torque() const{return torque(com());};
        flt moment_about(const Vec axis, const Vec loc) const;
        flt moment_about(const Vec axis) const{return moment_about(axis, com());};
        Vec angular_momentum(const Vec loc) const;
        Vec angular_momentum() const {return angular_momentum(com());};
        Matrix moment(const Vec loc) const;
        Matrix moment() const{return moment(com());};
        Vec omega(const Vec loc) const;
        Vec omega() const {return omega(com());};
        void add_omega(Vec w, Vec origin);
        void add_omega(Vec w){return add_omega(w, com());};
        inline void reset_L(){
            Vec c = com(), w = omega(c);
            if (w.squaredNorm() == 0) return;
            add_omega(-w, c);
        }
        #elif defined VEC2D
        flt torque(const Vec loc) const;
        flt moment(const Vec loc) const;
        flt moment() const{return moment(com());};
        flt angular_momentum(const Vec loc) const;
        flt angular_momentum() const {return angular_momentum(com());};
        flt omega(const Vec loc) const{return angular_momentum(loc) / moment(loc);};
        flt omega() const{return omega(com());};
        void add_omega(flt w, Vec origin);
        void add_omega(flt w){add_omega(w, com());}
        inline void reset_L(){
            Vec c = com();
            flt w = omega(c);
            if (w == 0) return;
            add_omega(-w, c);
        }
        #endif

        void add_velocity(Vec v);
        void reset_com_velocity(){add_velocity(-com_velocity());};
        void randomize_velocities(flt T);

        void reset_forces();
        virtual ~AtomGroup(){};
};

inline Atom& AtomIter::operator*() const{return g[i];};

class AtomVec : public virtual AtomGroup {
    private:
        Atom* atoms;
        uint sz;
        void zero() {
            for(uint i=0; i < sz; i++){
                atoms[i].x = Vec::Zero();
                atoms[i].v = Vec::Zero();
                atoms[i].f = Vec::Zero();
                atoms[i].a = Vec::Zero();
            }
        }
    public:
        AtomVec(vector<double> masses) : sz((uint) masses.size()){
            atoms = new Atom[sz];
            for(uint i=0; i < sz; i++) atoms[i].m = masses[i];
            zero();
        };
        AtomVec(uint N, flt mass) : sz(N){
            atoms = new Atom[sz];
            for(uint i=0; i < sz; i++) atoms[i].m = mass;
            zero();
        };
        AtomVec(AtomVec& other) : sz(other.size()){
            atoms = new Atom[sz];
            for(uint i=0; i < sz; i++) atoms[i] = other.atoms[i];
        };
        AtomVec& vec(){return *this;};
        inline Atom& operator[](cuint n){return atoms[n];};
        inline Atom& operator[](cuint n) const {return atoms[n];};
        inline AtomID get_id(cuint n) {
            if (n > sz) return AtomID(); return AtomID(atoms + n,n);};
        inline uint size() const {return sz;};

        ~AtomVec(){ delete [] atoms;};
};

class SubGroup : public AtomGroup {
    protected:
        sptr<AtomVec> atoms;
        vector<AtomID> ids;
    public:
        SubGroup(sptr<AtomVec> atoms) : atoms(atoms){};
        AtomVec &vec(){return *atoms;};
        inline Atom& operator[](cuint n){return *ids[n];};
        inline Atom& operator[](cuint n) const{return *ids[n];};
        inline Atom& get(cuint n){return *ids[n];};
        inline void add(AtomID a){
            std::vector<AtomID>::iterator it = std::find(ids.begin(), ids.end(), a);
            if(it != ids.end())
                throw std::invalid_argument("Cannot add AtomID to SubGroup: it already exists.");
            return ids.push_back(a);
        };
        inline AtomID get_id(cuint n) {return ids[n];};
        inline uint size() const {return (uint) ids.size();};
};

#endif