particle_group.h

Go to the documentation of this file.

#pragma once
 
#include <iomanip>
#include "Float.h"
#include "list.h"
 
namespace COMM {
 
 
    template <class Tparticle, class Tpcm>
    class ParticleGroup: public List<Tparticle>{
    private:
        typedef List<Tparticle> TList;
        bool origin_frame_flag; 
 
    public:
        Tpcm cm;  
 
 
        ParticleGroup(): TList(), origin_frame_flag(true) {};
  
 
        void clear() {
            TList::clear();
            origin_frame_flag = true;
        }
 
        /* Copy function will remove the local data and also copy the particle data or the link
         */
        ParticleGroup& operator = (const ParticleGroup& _particle_group) {
            TList::clear();
            *(List<Tparticle>*)this = *(List<Tparticle>*)&_particle_group;
            origin_frame_flag = _particle_group.origin_frame_flag;
            cm = _particle_group.cm;
 
            return *this;
        }
 
 
        ~ParticleGroup() {
            clear();
        }
 
 
        int getBackupDataSize() const {
            return 6*TList::num_;
        }
 
 
        int backupParticlePosVel(Float* _bk) {
            for (int i=0; i<TList::num_; i++) {
                const int k=6*i;
                Float* pos = TList::data_[i].getPos();
                Float* vel = TList::data_[i].getVel();
                _bk[k  ] = pos[0];
                _bk[k+1] = pos[1];
                _bk[k+2] = pos[2];
                _bk[k+3] = vel[0];
                _bk[k+4] = vel[1];
                _bk[k+5] = vel[2];
            }
            return 6*TList::num_;
        }
 
 
        int restoreParticlePosVel(Float* _bk) {
            for (int i=0; i<TList::num_; i++) {
                const int k=6*i;
                Float* pos = TList::data_[i].getPos();
                Float* vel = TList::data_[i].getVel();
                pos[0] = _bk[k  ];
                pos[1] = _bk[k+1];
                pos[2] = _bk[k+2];
                vel[0] = _bk[k+3];
                vel[1] = _bk[k+4];
                vel[2] = _bk[k+5];
            }
            return 6*TList::num_;
        }
 
 
        void printColumnTitle(std::ostream & _fout, const int _width=20) {
            _fout<<std::setw(_width)<<"N";
            for (int i=0; i<TList::num_; i++) TList::data_[i].printColumnTitle(_fout, _width);
        }
 
 
        void printColumn(std::ostream & _fout, const int _width=20){
            _fout<<std::setw(_width)<<TList::num_;
            for (int i=0; i<TList::num_; i++) TList::data_[i].printColumn(_fout, _width);
        }
 
 
        void writeBinary(FILE* _fout) {
            int num = TList::num_;
            fwrite(&num,sizeof(int), 1, _fout);
            for (int i=0; i<TList::num_; i++) TList::data_[i].writeBinary(_fout);
            fwrite(&origin_frame_flag, sizeof(bool), 1, _fout);
            cm.writeBinary(_fout);
        }
 
        //  @param [in] _fout: FILE IO for writing
        //*/
        //void writeAscii(FILE* _fout) const{
        //    fprintf(_fout, "%d ", TList::num_);
        //    for (int i=0; i<TList::num_; i++) TList::data_[i].writeAscii(_fout);
        //}
 
 
        void writeMemberAscii(std::ostream& _fout) const{
            _fout<<TList::num_<<" ";
            for (int i=0; i<TList::num_; i++) TList::data_[i].writeAscii(_fout);
        }
 
 
        void readBinary(FILE* _fin) {
            ASSERT(TList::mode_==ListMode::local);
            ASSERT(TList::num_==0);
            ASSERT(TList::nmax_==0);
            int n_new;
            int rn = fread(&n_new, sizeof(int),1, _fin);
            if(rn<1) {
                std::cerr<<"Error: cannot read particle number!\n";
                abort();
            }
            if(n_new<=0) {
                std::cerr<<"Error: reading particle number "<<n_new<<"<=0!\n";
                abort();
            }
            TList::reserveMem(n_new);
            for (int i=0; i<n_new; i++) TList::data_[i].readBinary(_fin);
            TList::num_ = n_new;
            rn = fread(&origin_frame_flag, sizeof(bool), 1, _fin);
            if(rn<1) {
                std::cerr<<"Error: cannot read origin_frame_flag!\n";
                abort();
            }
            cm.readBinary(_fin);
        }
 
 
        void readMemberAscii(std::istream& _fin) {
            ASSERT(TList::mode_==ListMode::local);
            ASSERT(TList::num_==0);
            ASSERT(TList::nmax_==0);
            int n_new;
            _fin>>n_new;
            ASSERT(!_fin.eof());
            if(n_new<=0) {
                std::cerr<<"Error: reading particle number "<<n_new<<"<=0!\n";
                abort();
            }
            TList::reserveMem(n_new);
            for (int i=0; i<n_new; i++) {
                TList::data_[i].readAscii(_fin);
                ASSERT(!_fin.eof());
            }
            TList::num_ = n_new;
        }
 
 
        void shiftToCenterOfMassFrame() {
            if (origin_frame_flag) {
                const Float *rc = cm.getPos();
                const Float *vc = cm.getVel();
                for (int i=0;i<TList::num_;i++) {
                    Float *ri = TList::data_[i].getPos();
                    Float *vi = TList::data_[i].getVel();
                    ri[0] -= rc[0];
                    ri[1] -= rc[1];
                    ri[2] -= rc[2];
                    vi[0] -= vc[0];
                    vi[1] -= vc[1];
                    vi[2] -= vc[2];
                }
                origin_frame_flag = false;
            }
            else {
                std::cerr<<"Warning: particles are already in the center-of-mass frame!\n";
            }      
        }
        
 
        void shiftToOriginFrame() {
            if (origin_frame_flag) {
                std::cerr<<"Warning: particles are already in original frame!\n";
            }
            else {
                const Float *rc = cm.getPos();
                const Float *vc = cm.getVel();
                for (int i=0;i<TList::num_;i++) {
                    Float *ri = TList::data_[i].getPos();
                    Float *vi = TList::data_[i].getVel();
                    ri[0] += rc[0];
                    ri[1] += rc[1];
                    ri[2] += rc[2];
                    vi[0] += vc[0];
                    vi[1] += vc[1];
                    vi[2] += vc[2];
                }
                origin_frame_flag = true;
            }
        }
 
        void calcCenterOfMass() {
            Float *rc = cm.getPos();
            Float *vc = cm.getVel();
            rc[0] = rc[1] = rc[2] = 0.0;
            vc[0] = vc[1] = vc[2] = 0.0;
            cm.mass = 0.0;
            for (int i=0;i<TList::num_;i++) {
                const Float *ri = TList::data_[i].getPos();
                const Float *vi = TList::data_[i].getVel();
                const Float mi = TList::data_[i].mass;
                rc[0] += ri[0] * mi;
                rc[1] += ri[1] * mi;
                rc[2] += ri[2] * mi;
 
                vc[0] += vi[0] * mi;
                vc[1] += vi[1] * mi;
                vc[2] += vi[2] * mi;
 
                cm.mass += mi;
            }
            rc[0] /= cm.mass; 
            rc[1] /= cm.mass; 
            rc[2] /= cm.mass; 
            vc[0] /= cm.mass; 
            vc[1] /= cm.mass; 
            vc[2] /= cm.mass;
        }
        
        bool isOriginFrame() const {
            return origin_frame_flag;
        }
    };
 
}