StatisticsVector.h

Go to the documentation of this file.

#ifndef FSTAT_H
#define FSTAT_H
#include "MD.h"
#include "Matrix.h"
#include <string.h>
enum StatType {O, X, Y, Z, XY, XZ, YZ, XYZ};
#include "Polynomial.h"
#include "StatisticsPoint.h"

template <StatType T>
class StatisticsVector : public virtual MD
{
public:
        void constructor();
        void constructor(string name);
        
        StatisticsVector() {constructor();}
        StatisticsVector(string name) {constructor(name);}

        StatisticsVector(StatisticsVector& other);
                
        StatisticsVector(unsigned int argc, char *argv[]);
        void readStatArguments(unsigned int argc, char *argv[]);

        string print();

        void set_statType();

        void print_help();

        //set functions should be called before executing the statistical routine 
        void set_nx(int new_) { nx = new_;};
        void set_hx(Mdouble hx) { set_nx( ceil((get_xmaxStat()-get_xminStat())/hx) ); };
        void set_ny(int new_) { ny = new_;};
        void set_hy(Mdouble hy) { set_ny( ceil((get_ymaxStat()-get_yminStat())/hy) ); };
        void set_nz(int new_) { nz = new_; if (get_dim()<3) cerr << "Warning in set_nz: dimension less than 3"<<endl; };
        void set_hz(Mdouble hz) { set_nz( ceil((get_zmaxStat()-get_zminStat())/hz) ); };
        void set_n(int n) { set_nx(n); set_ny(n); set_nz(n); };
        void set_h(Mdouble h) { set_hx(h); set_hy(h); set_hz(h);};
        int get_nx() {return nx;};
        int get_ny() {return ny;};
        int get_nz() {return nz;};
        void set_tminStat(Mdouble t) {tminStat=t;};
        void set_tmaxStat(Mdouble t) {tmaxStat=t;};
        void set_tintStat(Mdouble t) {tintStat=t;};
        Mdouble get_tminStat() {return tminStat;};
        Mdouble get_tmaxStat() {if(isnan(tmaxStat)) return get_tmax(); else return tmaxStat;};
        Mdouble get_tintStat() {return tintStat;};
        bool check_current_time_for_statistics(){return (get_t()>get_tminStat()&&get_t()<=get_tmaxStat()+get_dt());};
        
        void set_CG_type(const char* new_);
        void set_CG_type(CG new_);
        CG get_CG_type() {return CG_type;};
        
        void set_n(int nx_, int ny_, int nz_) {nx=nx_; ny=ny_; nz=nz_;};
        void get_n(int& nx_, int& ny_, int& nz_) {nx_=nx; ny_=ny; nz_=nz;};

        void set_w(Mdouble w) {set_w2(sqr(w));};

        void set_w2(Mdouble new_);

        Mdouble get_w() {return sqrt(w2);};
        Mdouble get_w2() {return w2;};
        Mdouble get_cutoff() {return cutoff;};
        Mdouble get_cutoff2() {return sqr(cutoff);};
        
        //~ bool get_boundedDomain() {return boundedDomain;};
        //~ void set_boundedDomain(bool new_) {boundedDomain=new_;};

        string print_CG();

        StatisticsPoint<T> average(vector<StatisticsPoint<T> > &P);

        void reset_statistics();

        void statistics_from_fstat_and_data();
        
        void set_doTimeAverage(bool new_) {doTimeAverage = new_;};
        bool get_doTimeAverage() {return doTimeAverage;};
        
        void set_StressTypeForFixedParticles(int new_) {StressTypeForFixedParticles = new_;};
        int get_StressTypeForFixedParticles() {return StressTypeForFixedParticles;};
        
        // to keep compatible with older versions
        void set_infiniteStressForFixedParticles(bool new_) {StressTypeForFixedParticles = 2+new_;};
        
        void set_mirrorAtDomainBoundary(Mdouble new_) {mirrorAtDomainBoundary = new_;};
        Mdouble get_mirrorAtDomainBoundary() {return mirrorAtDomainBoundary;};
        
        void set_doVariance(bool new_) {doVariance = new_;};
        bool get_doVariance() {return doVariance;};

        void set_doGradient(bool new_) {doGradient = new_;};
        bool get_doGradient() {return doGradient;};
        
        void set_superexact(bool new_) {superexact = new_;};
        bool get_superexact() {return superexact;};
        
        void set_ignoreFixedParticles(bool new_) {ignoreFixedParticles = new_;};
        bool get_ignoreFixedParticles() {return ignoreFixedParticles;};
        
        void verbose() {verbosity = 2;};
        void set_verbosity(int new_) {verbosity = new_;};
        int get_verbosity() {return verbosity;};
        
        void set_walls(bool new_) {walls = new_;};
        bool get_walls() {return walls;};

        void set_periodicWalls(bool new_) {periodicWalls = new_;};
        bool get_periodicWalls() {return periodicWalls;};

        void set_w_over_rmax(Mdouble new_) {w_over_rmax = new_;};
        Mdouble get_w_over_rmax() {return w_over_rmax;};

        void set_Positions();

        bool read_next_from_data_file (unsigned int format);

        void gather_force_statistics_from_fstat_and_data();
        void evaluate_force_statistics(int wp=0);
        void evaluate_wall_force_statistics(Vec3D P, int wp=0);
        void jump_fstat();
        
        void initialize_statistics();
        
        void output_statistics();
        void gather_statistics_collision(int index1,int index2, Vec3D Contact,  Mdouble delta, Mdouble ctheta, Mdouble fdotn, Mdouble fdott, Vec3D P1_P2_normal_, Vec3D P1_P2_tangential);
        void process_statistics(bool usethese);
        void finish_statistics();

        void write_statistics();
        void write_time_average_statistics();   
        
        void evaluate_particle_statistics(vector<Particle*>::iterator P, int wp=0);
        
        vector<StatisticsPoint<T> > get_Points(){return Points;};
        
        Mdouble get_xminStat() {if(isnan(xminStat)) return get_xmin(); else return xminStat;};
        Mdouble get_yminStat() {if(isnan(yminStat)) return get_ymin(); else return yminStat;};
        Mdouble get_zminStat() {if(isnan(zminStat)) return get_zmin(); else return zminStat;};
        Mdouble get_xmaxStat() {if(isnan(xmaxStat)) return get_xmax(); else return xmaxStat;};
        Mdouble get_ymaxStat() {if(isnan(ymaxStat)) return get_ymax(); else return ymaxStat;};
        Mdouble get_zmaxStat() {if(isnan(zmaxStat)) return get_zmax(); else return zmaxStat;};
        void set_xminStat(Mdouble xminStat_){xminStat=xminStat_;};
        void set_yminStat(Mdouble yminStat_){yminStat=yminStat_;};
        void set_zminStat(Mdouble zminStat_){zminStat=zminStat_;};
        void set_xmaxStat(Mdouble xmaxStat_){xmaxStat=xmaxStat_;};
        void set_ymaxStat(Mdouble ymaxStat_){ymaxStat=ymaxStat_;};
        void set_zmaxStat(Mdouble zmaxStat_){zmaxStat=zmaxStat_;};
        int get_nTimeAverage() {return nTimeAverage;};

        Mdouble setInfinitelyLongDistance();

        void set_Polynomial(vector<Mdouble> new_coefficients, unsigned int new_dim) {   
                Polynomial.set_polynomial(new_coefficients, new_dim);
        }

        void set_Polynomial(Mdouble* new_coefficients, unsigned int num_coeff, unsigned int new_dim) {  
                Polynomial.set_polynomial(new_coefficients, num_coeff, new_dim);
        }

        void set_PolynomialName(const char* new_name) { 
                Polynomial.set_name(new_name);
        }

        string get_PolynomialName() {   
                return Polynomial.get_name();
        }

        void set_rmin(Mdouble new_) {rmin=new_;}
        void set_rmax(Mdouble new_) {rmax=new_;}
        void set_hmax(Mdouble new_) {hmax=new_;}

        Mdouble evaluatePolynomial(Mdouble r) {
                return Polynomial.evaluate(r);
        }

        Mdouble evaluateIntegral(Mdouble n1, Mdouble n2, Mdouble t) {
                return Polynomial.evaluateIntegral(n1,n2,t);
        }


//Member Variables
private:

        //General Variables
        StatType statType;
        int nx;
        int ny;
        int nz; 
        Mdouble xminStat, xmaxStat, yminStat, ymaxStat, zminStat, zmaxStat;
        int nxMirrored, nyMirrored, nzMirrored; 
        
        //Storage of points and gradients
        vector<StatisticsPoint<T> > Points;
        vector<StatisticsPoint<T> > dx; 
        vector<StatisticsPoint<T> > dy; 
        vector<StatisticsPoint<T> > dz; 

        //For time averaging
        vector<StatisticsPoint<T> > timeAverage; 
        vector<StatisticsPoint<T> > timeVariance; 
        vector<StatisticsPoint<T> > dxTimeAverage; 
        vector<StatisticsPoint<T> > dyTimeAverage; 
        vector<StatisticsPoint<T> > dzTimeAverage; 
        bool doTimeAverage;
        int nTimeAverageReset; 
        bool doVariance;
        bool doGradient;
        int nTimeAverage; 
        
        //Coarse graining variables
        CG CG_type; 
        NORMALIZED_POLYNOMIAL<T> Polynomial;
        // ///<if true, then the course-graining function will be cut at the domain boundaries and resized to satisfy int(W) = 1
        // bool boundedDomain; 
        Mdouble w2; 
        Mdouble cutoff, cutoff2; 
        Mdouble w_over_rmax; 
        
        //Options that can be set before evaluation

        Mdouble tminStat; 
        Mdouble tmaxStat; 
        Mdouble tintStat; 
        Mdouble indSpecies; 
        Mdouble rmin; 
        Mdouble rmax; 
        Mdouble hmax; 
        bool walls; 
        bool periodicWalls; 
        bool ignoreFixedParticles; 
        int StressTypeForFixedParticles;

        int verbosity; //<Determines how much is outputted to the terminal 
        int format; //< format of the data input file
        Mdouble mirrorAtDomainBoundary; //< 0: Statistics near the wall are equal to statistics away from the wall; 1: Statistics are mirrored at the domain boundaries; up to a maximum depth specified by this number

        bool isMDCLR;
        bool superexact; 
        
        


        //Variables communicate values between member functions #evaluate_force_statistics and #gather_statistics_collision)used to communicate values between member functions evaluate_force_statistics and gather_force_statistics   
        Vec3D P1; 
        Vec3D P2; 
        Vec3D P1_P2_normal; 
        Mdouble P1_P2_distance; 
        Matrix3D P1_P2_NormalStress; 
        Matrix3D P1_P2_TangentialStress; 
        Vec3D P1_P2_NormalTraction; 
        Vec3D P1_P2_TangentialTraction; 
        MatrixSymmetric3D P1_P2_Fabric; 
        Vec3D P1_P2_CollisionalHeatFlux; 
        Mdouble P1_P2_Dissipation; 
        Mdouble P1_P2_Potential; 

};

#include "StatisticsPoint.hcc"
#include "StatisticsVector.hcc"

#endif