StatisticsVector.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 FSTAT_H
#define FSTAT_H
#include "MD.h"
#include "MatrixSymmetric.h"
#include <string.h>
#include <fstream>
#include <math.h>
enum StatType {O, X, Y, Z, XY, XZ, YZ, XYZ, RAZ, RZ, AZ, RA, R, A};
#include "Polynomial.h"
#include "StatisticsPoint.h"

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

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

    std::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) std::cerr << "Warning in set_nz: dimension less than 3"<<std::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(std::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_;};

    std::string print_CG();

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

    virtual void reset_statistics();

    void statistics_from_fstat_and_data();
    void statistics_from_p3();
    void jump_p3c();
    
    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 gather_force_statistics_from_p3c(int version);
    void gather_force_statistics_from_p3w(int version, std::vector<int>& index);
    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(std::vector<BaseParticle*>::iterator P, int wp=0);
    
    std::vector<StatisticsPoint<T> > get_Points(){return Points;};
    
    Mdouble get_xminStat() {if(std::isnan(xminStat)) return get_xmin(); else return xminStat;};
    Mdouble get_yminStat() {if(std::isnan(yminStat)) return get_ymin(); else return yminStat;};
    Mdouble get_zminStat() {if(std::isnan(zminStat)) return get_zmin(); else return zminStat;};
    Mdouble get_xmaxStat() {if(std::isnan(xmaxStat)) return get_xmax(); else return xmaxStat;};
    Mdouble get_ymaxStat() {if(std::isnan(ymaxStat)) return get_ymax(); else return ymaxStat;};
    Mdouble get_zmaxStat() {if(std::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(std::vector<Mdouble> new_coefficients, unsigned int new_dim) {  
        CGPolynomial.set_polynomial(new_coefficients, new_dim);
    }

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

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

    std::string get_PolynomialName() {  
        return CGPolynomial.get_name();
    }

    void set_doublePoints(bool new_) {doublePoints=new_;}

    bool get_doublePoints() {return doublePoints;}

    void set_TimeAverageReset(int new_) {nTimeAverageReset=new_;}

    bool get_TimeAverageReset() {return nTimeAverageReset;}

    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 CGPolynomial.evaluate(r);
    }

    Mdouble evaluatePolynomialGradient(Mdouble r) {
        return CGPolynomial.evaluateGradient(r);
    }

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


//Member Variables
protected:

    //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
    std::vector<StatisticsPoint<T> > Points;
    std::vector<StatisticsPoint<T> > dx; 
    std::vector<StatisticsPoint<T> > dy; 
    std::vector<StatisticsPoint<T> > dz; 

    //For time averaging
    std::vector<StatisticsPoint<T> > timeAverage; 
    std::vector<StatisticsPoint<T> > timeVariance; 
    std::vector<StatisticsPoint<T> > dxTimeAverage; 
    std::vector<StatisticsPoint<T> > dyTimeAverage; 
    std::vector<StatisticsPoint<T> > dzTimeAverage; 
    bool doTimeAverage;
    int nTimeAverageReset; 
    bool doVariance;
    bool doGradient;
    int nTimeAverage; 
    
    //Coarse graining variables
    CG CG_type; 
    NORMALIZED_POLYNOMIAL<T> CGPolynomial;
    // ///<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; 
    
    //uses close points to allow calculation of gradients
    bool doublePoints; 
    
    bool satisfiesInclusionCriteria(BaseParticle* P);


    //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_ContactCoupleStress; 
    Vec3D P1_P2_Contact; 
    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; 

    bool loadVelocityProfile (const char* filename);
    Vec3D getVelocityProfile (Vec3D Position);

    std::vector<Vec3D> VelocityProfile;
    Vec3D VelocityProfile_Min;
    Vec3D VelocityProfile_D;
    
    std::fstream p3p_file;
    std::fstream p3c_file;
    std::fstream p3w_file;
    bool read_next_from_p3p_file ();
    void auto_setdim (); 
    


};

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

#endif