StatisticsVector.hcc

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.

std::ostream& operator<<(std::ostream& os, const StatType S) {
    if (S==O) os << "O";
    else if (S==X) os << "X";
    else if (S==Y) os << "Y";
    else if (S==Z) os << "Z";
    else if (S==XY) os << "XY";
    else if (S==XZ) os << "XZ";
    else if (S==YZ) os << "YZ";
    else if (S==XYZ) os << "XYZ";
    else if (S==RAZ) os << "RAZ";
    else if (S==RA) os << "RA";
    else if (S==AZ) os << "AZ";
    else if (S==RZ) os << "RZ";
    else if (S==R) os << "R";
    else if (S==A) os << "A";
    return os;
}

template <StatType T>
void StatisticsVector<T>::set_w2(Mdouble new_) {
    // set w2
    if (new_>0.0) w2 = new_;
    else w2 = sqr(w_over_rmax*getLargestParticle()->get_Radius());
    //(2.0*get_xmax()/nx);
    
    // set cutoff radius
    if (get_CG_type()==HeavisideSphere||get_CG_type()==Polynomial) {
        cutoff = sqrt(w2);
    } else if (get_CG_type()==Gaussian) {
        if (get_superexact()) cutoff = 5.0*sqrt(w2);
        else cutoff = 3.0*sqrt(w2);
    } else { std::cerr << "error in CG_function" << std::endl; exit(-1); }
}

template <StatType T>
void StatisticsVector<T>::reset_statistics() {
    for (unsigned int i=0; i<Points.size(); i++) Points[i].set_zero(); 
    if (get_doGradient()) for (unsigned int i=0; i<Points.size(); i++) {
        dx[i].set_zero(); dy[i].set_zero(); dz[i].set_zero(); 
    }
    for (std::vector<BaseParticle*>::iterator it = getParticleHandler().begin(); it!=getParticleHandler().end(); ++it) {
        (*it)->set_PreviousPosition((*it)->get_Position());
    }
}

template <StatType T>
void StatisticsVector<T>::constructor() { 
    //stores the template parameter Stattype in a variable
    set_statType();
    
    // sets connection between StatisticsPoint and StatisticsVector
    StatisticsPoint<T>::set_gb(this);
    // set nx, ny,nz
    nx = ny = nz = 1;
    nxMirrored=nyMirrored=nzMirrored=0;

    // set default CG variables 
    CG_type = Gaussian;
    // in get_statistics_..., w2 is set to w_over_rmax*rmax if w2==0
    w_over_rmax = 1;
    
    // bounded domain
    //~ boundedDomain = false;
    
    // set default values; variables will only be used if the default 
    // is overwritten; thus the default is a nonsense value 
    w2 = 0.0;
    set_tintStat(0);
    set_tminStat(-1e20);
    set_tmaxStat(NAN);
    rmin = 0;
    hmax = 1e20;
    rmax = 1e20;
    indSpecies = -1; //all species

    
    //set default options - Not used it physially don't make sense i.e if ?max<?min
    xminStat = NAN;
    xmaxStat = NAN;
    yminStat = NAN;
    ymaxStat = NAN;
    zminStat = NAN;
    zmaxStat = NAN; 
    

    // set default options
    ignoreFixedParticles = true; //this should be true if you want statistics only of the flow
    StressTypeForFixedParticles = 1;
    verbosity = 1;
    walls = true;
    periodicWalls = true;
    format = 0;
    mirrorAtDomainBoundary = false;
    set_superexact(false);
    VelocityProfile.resize(0); //equivalent to do not use velocityprofile

    // time averaging
    set_doTimeAverage(true);
    nTimeAverage = 0;
    nTimeAverageReset = -1;
    //calculate variance
    set_doVariance(false);
    //calculate gradient
    set_doGradient(false);
    doublePoints=false;

    // additional stuff
    set_options_stat(1);
}

template <StatType T>
StatisticsVector<T>::StatisticsVector(StatisticsVector& other): MD()
{
    //assumes that we want the statistical method copied, but resets the time averaging 
    constructor();
    //~ StatType statType; //this has to be constant, right?
    nx = other.nx;
    ny = other.ny;
    nz = other.nz;
    xminStat = other.xminStat;
    yminStat = other.yminStat;
    zminStat = other.zminStat;
    xmaxStat = other.xmaxStat;
    ymaxStat = other.ymaxStat;
    zmaxStat = other.zmaxStat;
    nxMirrored = other.nxMirrored;
    nyMirrored = other.nyMirrored;
    nzMirrored = other.nzMirrored;
    Points = other.Points;
    dx = other.dx;
    dy = other.dy;
    dz = other.dz;

    timeAverage = other.timeAverage;
    timeVariance = other.timeVariance;
    dxTimeAverage = other.dxTimeAverage;
    dyTimeAverage = other.dyTimeAverage;
    dzTimeAverage = other.dzTimeAverage;
    doTimeAverage = other.doTimeAverage;
    nTimeAverage = 0;
    for (unsigned int i=0; i<timeAverage.size(); i++) timeAverage[i].set_zero();
    for (unsigned int i=0; i<timeVariance.size(); i++) timeVariance[i].set_zero();
    for (unsigned int i=0; i<dxTimeAverage.size(); i++) dxTimeAverage[i].set_zero();
    for (unsigned int i=0; i<dyTimeAverage.size(); i++) dyTimeAverage[i].set_zero();
    for (unsigned int i=0; i<dzTimeAverage.size(); i++) dzTimeAverage[i].set_zero();
    doVariance = other.doVariance;
    doGradient = other.doGradient;
    CG_type = other.CG_type;
    CGPolynomial = other.CGPolynomial;
    w2 = other.w2;
    cutoff = other.cutoff;
    cutoff2 = other.cutoff2;
    w_over_rmax = other.w_over_rmax;
    tminStat = other.tminStat;
    tmaxStat = other.tmaxStat;
    tintStat = other.tintStat;
    rmin = other.rmin;
    rmax = other.rmax;
    hmax = other.hmax;
    indSpecies = other.indSpecies;
    walls = other.walls;
    periodicWalls = other.periodicWalls;
    ignoreFixedParticles = other.ignoreFixedParticles;
    StressTypeForFixedParticles = other.StressTypeForFixedParticles;
    verbosity = other.verbosity;
    format = other.format;
    mirrorAtDomainBoundary = other.mirrorAtDomainBoundary;
    isMDCLR = other.isMDCLR;
    superexact = other.superexact;
}



template <StatType T>
void StatisticsVector<T>::constructor(std::string name) { 
    //call default constructor
    constructor();
    
    // set name
    set_name(name.c_str()); 

    if (!strcmp(get_name().c_str(),"c3d")) {
        //write here what to do with Stefan's data 
        std::cout << "MDCLR data" << std::endl;
        exit(-1);
    } else {
        //write here what to do with Mercury data
        isMDCLR = false;
        load_restart_data();
        set_append(false);
        //~ set_tmaxStat(t+dt); //note:doesn't work if restart data is not the latest
    }
}

//Constructor that accepts arguments from the command line
template <StatType T>
StatisticsVector<T>::StatisticsVector(unsigned int argc, char *argv[])
{
    //check if filename is given
    if (argc<2) {
        std::cerr<<"fstatistics.exe filename [-options]"<< std::endl
        << "Type -help for more information" << std::endl;
        exit(-1);
    } 
    
    //print help if required
    if (!strcmp(argv[1],"-help")) {
        print_help(); 
    }

    //create StatisticsVector
    std::string filename(argv[1]);
    constructor(filename);

    //check for options (standardly '-option argument')
    readStatArguments(argc, argv);
}

template <StatType T>
bool StatisticsVector<T>::loadVelocityProfile (const char* filename) 
{
    //This opens the file the data will be recalled from
    std::fstream is;
    is.open(filename, std::fstream::in);
    if (!is.is_open()||is.bad()) {
        std::cout << "Loading data file " << filename << " failed" << std::endl;
        return false;
    } 

    //Read the file; we assume that the stat file we read in has the same stattype and domain size as the problem.
    std::string line_string;
    getline(is,line_string);
    std::string dummy;
    is >> dummy >> dummy;
    is >> dummy >> dummy;
    double xmin, xmax, ymin, ymax, zmin, zmax;
    is >> dummy >> xmin >> xmax >> ymin >> ymax >> zmin >> zmax;
    VelocityProfile_Min = Vec3D(xmin,ymin,zmin);
    double nx, ny, nz;
    is >> dummy >> nx >> ny >> nz;
    double dx = (xmax-xmin)/(nx-1);
    double dy = (ymax-ymin)/(ny-1);
    double dz = (zmax-zmin)/(nz-1);
    VelocityProfile_D = Vec3D(dx,dy,dz);
    std::string statType;
    is >> dummy >> statType;
    getline(is,line_string);
    getline(is,line_string);
    
    set_nx(nx);
    set_ny(ny);
    set_nz(nz);
    set_xminStat(xmin);
    set_yminStat(ymin);
    set_zminStat(zmin);
    set_xmaxStat(xmax);
    set_ymaxStat(ymax);
    set_zmaxStat(zmax);
    
    VelocityProfile.resize(nx*ny*nz);
    Vec3D M;
    double rho;
    int n=0;
    for (int i=0; i<nx; i++)
    for (int j=0; j<ny; j++)
    for (int k=0; k<nz; k++) {
        is >> dummy >> dummy >> dummy >> dummy >> rho >> M.X >> M.Y >> M.Z;
        getline(is,line_string);
        VelocityProfile[n]=M/rho;
        n++;
    }
    
    //Close the file
    is.close();
    std::cout << "Loaded data file " << filename << std::endl;
    return true;
}

template <StatType T>
Vec3D StatisticsVector<T>::getVelocityProfile (Vec3D Position) 
{
    double i,j,k;
    double l,m,n;
    l = modf ((double)(Position.X-VelocityProfile_Min.X)/VelocityProfile_D.X, &i);  
    m = modf ((double)(Position.Y-VelocityProfile_Min.Y)/VelocityProfile_D.Y, &j);  
    n = modf ((double)(Position.Z-VelocityProfile_Min.Z)/VelocityProfile_D.Z, &k);
    if (nx<1) i=0;
    if (ny<1) j=0;
    if (ny<1) k=0;
    int ix = (i*ny+j)*nz+k;
    Vec3D V;
    if (nx<=1) V.X = VelocityProfile[ix].X;
    else V.X=VelocityProfile[ix].X*(1-l)+VelocityProfile[ix+ny*nz].X*l;
    if (ny<=1) V.Y = VelocityProfile[ix].Y;
    else V.Y=VelocityProfile[ix].Y*(1-m)+VelocityProfile[ix+nz].Y*m;
    if (nz<=1) V.Z = VelocityProfile[ix].Z;
    else V.Z=VelocityProfile[ix].Z*(1-n)+VelocityProfile[ix+1].Z*n;
    return V;
}


//Constructor that accepts arguments from the command line
template <StatType T>
void StatisticsVector<T>::readStatArguments(unsigned int argc, char *argv[])
{
    //check for options (standardly '-option argument')
    for (unsigned int i = 2; i<argc; i+=2) {
        std::cout << "interpreting input argument " << argv[i];
        if (i+1<argc) std::cout << " " << argv[i+1];
        std::cout << std::endl;
        
        if (!strcmp(argv[i],"-CGtype")) {
            //CG_type_todo
            if (!strcmp(argv[i+1],"HeavisideSphere")) {
                set_CG_type(HeavisideSphere);
            } else if (!strcmp(argv[i+1],"Gaussian")) {
                set_CG_type(Gaussian);
            } else {
                set_CG_type(argv[i+1]);
            } 
        } else if (!strcmp(argv[i],"-w")) {
            double old = get_w();
            set_w (atof(argv[i+1]));
            std::cout << " w changed from " << old << " to " << get_w() << std::endl;
        } else if (!strcmp(argv[i],"-help")) {
            print_help(); 
            i--; //requires no argument
        } else if (!strcmp(argv[i],"-velocityprofile")) {
            loadVelocityProfile(argv[i+1]);
        } else if (!strcmp(argv[i],"-h")) {
            set_h (atof(argv[i+1]));
        } else if (!strcmp(argv[i],"-hx")) {
            set_hx(atof(argv[i+1]));
        } else if (!strcmp(argv[i],"-hy")) {
            set_hy(atof(argv[i+1]));
        } else if (!strcmp(argv[i],"-hz")) {
            set_hz(atof(argv[i+1]));
        } else if (!strcmp(argv[i],"-n")) {
            set_n(atoi(argv[i+1]));
        } else if (!strcmp(argv[i],"-nx")) {
            set_nx(atoi(argv[i+1]));
        } else if (!strcmp(argv[i],"-ny")) {
            set_ny(atoi(argv[i+1]));
        } else if (!strcmp(argv[i],"-nz")) {
            set_nz(atoi(argv[i+1]));
        } else if (!strcmp(argv[i],"-x")) {
            set_xminStat(atof(argv[i+1]));
            set_xmaxStat(atof(argv[i+2]));
            i++; //requires two arguments
        } else if (!strcmp(argv[i],"-y")) {
            set_yminStat(atof(argv[i+1]));
            set_ymaxStat(atof(argv[i+2]));
            i++; //requires two arguments
        } else if (!strcmp(argv[i],"-z")) {
            set_zminStat(atof(argv[i+1]));
            set_zmaxStat(atof(argv[i+2]));
            i++; //requires two arguments
        } else if (!strcmp(argv[i],"-indSpecies")) {
            indSpecies = atoi(argv[i+1]);
        } else if (!strcmp(argv[i],"-rmin")) {
            rmin = atof(argv[i+1]);
        } else if (!strcmp(argv[i],"-rmax")) {
            rmax = atof(argv[i+1]);
        } else if (!strcmp(argv[i],"-hmax")) {
            hmax = atof(argv[i+1]);
        } else if (!strcmp(argv[i],"-periodicwalls")) {
            set_periodicWalls(atoi(argv[i+1]));
        } else if (!strcmp(argv[i],"-walls")) {
            set_walls(atoi(argv[i+1]));
        } else if (!strcmp(argv[i],"-verbosity")) {
            set_verbosity(atoi(argv[i+1]));
        } else if (!strcmp(argv[i],"-verbose")) {
            verbose();
            i--; //requires no argument
        } else if (!strcmp(argv[i],"-format")) {
            format = atoi(argv[i+1]);
        } else if (!strcmp(argv[i],"-doublePoints")) {
            if (argc<=i+1||argv[i+1][0]=='-') { 
                set_doublePoints(true); 
                i--; //no argument
            } else set_doublePoints(atoi(argv[i+1]));
            std::cout << "set doublePoints=" << get_doublePoints() << std::endl;
        } else if (!strcmp(argv[i],"-w_over_rmax")) {
            set_w_over_rmax(atof(argv[i+1]));
        } else if (!strcmp(argv[i],"-tmin")) {
            set_tminStat(atof(argv[i+1]));
        } else if (!strcmp(argv[i],"-tmax")) {
            set_tmaxStat(atof(argv[i+1]));
        } else if (!strcmp(argv[i],"-tint")) {
            set_tintStat(atof(argv[i+1]));
        //~ } else if (!strcmp(argv[i],"-timesteps")) {
            //~ set_tintStat(atof(argv[i+1])*get_dt());
        } else if (!strcmp(argv[i],"-stepsize")) {
            int step_size = atoi(argv[i+1]);
            set_step_size(step_size);
        } else if (!strcmp(argv[i],"-loaddatafile")) {
            load_from_data_file(argv[i+1]);
            set_tminStat(get_t());
            set_tmaxStat(get_tmax());
        } else if (!strcmp(argv[i],"-statfilename") | !strcmp(argv[i],"-o")) {
            stat_filename.str(argv[i+1]);
            // set_name(argv[i+1]);
        } else if (!strcmp(argv[i],"-timeaverage")) {
            set_doTimeAverage(atoi(argv[i+1]));
        } else if (!strcmp(argv[i],"-timeaveragereset")) {
            nTimeAverageReset = atoi(argv[i+1]);
        } else if (!strcmp(argv[i],"-gradient")) {
            // use default argument if no argument is given
            if (argc<=i+1||argv[i+1][0]=='-') { 
                set_doGradient(true); 
                i--; //no argument
            } else set_doGradient(atoi(argv[i+1])); 
        } else if (!strcmp(argv[i],"-timevariance")) {
            // use default argument if no argument is given
            if (argc<=i+1||argv[i+1][0]=='-') { 
                set_doVariance(true); 
                i--; //no argument
            } else set_doVariance(atoi(argv[i+1])); 
        } else if (!strcmp(argv[i],"-superexact")) {
            // use default argument if no argument is given
            if (argc<=i+1||argv[i+1][0]=='-') { 
                set_superexact(true); 
                i--; //no argument
            } else set_superexact(atoi(argv[i+1])); 
        } else if (!strcmp(argv[i],"-ignoreFixedParticles")) {
            // use default argument if no argument is given
            if (argc<=i+1||argv[i+1][0]=='-') { 
                set_ignoreFixedParticles(true); 
                i--; //no argument
            } else set_ignoreFixedParticles(atoi(argv[i+1])); 
        } else if (!strcmp(argv[i],"-StressTypeForFixedParticles")) {
            StressTypeForFixedParticles = atoi(argv[i+1]);
        } else if (!strcmp(argv[i],"-mirrorAtDomainBoundary")) {
            set_mirrorAtDomainBoundary(atof(argv[i+1]));
        } else if (!strcmp(argv[i],"-stattype")) {
            //this was already used in Statistics()
        } else if (readNextArgument(i, argc, argv)) {
            //~ std::cout << " (read as MD argument, not statistics)" << std::endl; 
        } else {
            std::cerr << "Error: option unknown: " << argv[i] << std::endl;
            exit(-1);
        }
    }
}

//Constructor that accepts arguments from the command line
template <StatType T>
void StatisticsVector<T>::print_help()
{
    std::cout << "fstatistics.exe filename [-options]: " << std::endl
              << "This program evaluates statistical data from .fstat, .data, and. restart files and writes it into a .stat file" << std::endl << std::endl;
    std::cout << "OPTIONS:" << std::endl << std::endl;
    std::cout << "-stattype [X,Y,Z,XY,XZ,YZ,XYZ,RAZ,RA,AZ,RZ,R,A]: " << std::endl
             << "Used to obtain statistics averaged in all coordinate directions but the ones specified; f.e. -stattype Z yields depth-profiles. Default is XYZ." << std::endl << std::endl
         << "-CGtype [HeavisideSphere,Gaussian]: " << std::endl
             << "Averaging function; default: Gaussian" << std::endl << std::endl
         << "-w,-w_over_rmax [Mdouble]: " << std::endl
             << "Averaging width, absolute or in multiples of the radius of the largest particle in the restart file; default: w_over_rmax=1" << std::endl << std::endl
         << "-n,nx,ny,nz [uint]: " << std::endl
             << "Specifies the amount of grid points in coordinate directions for which statistics are evaluated; use -n to set all 3 directions at once; is set to one in averaged directions (see stattype). Default n=1" << std::endl << std::endl
         << "-h,hx,hy,hz [Mdouble]: " << std::endl
             << "Defines the amount of grid points by setting the mesh size " << std::endl << std::endl
         << "-x,y,z [Mdouble Mdouble]: " << std::endl
             << "Specifies the domain of the grid; f.e. for -x 0 1 all grid points will have x-values within [0,1]. Default -x xmin xmax (from restart file)" << std::endl << std::endl
         << "-indSpecies [int]: " << std::endl
         << "-rmin,rmax [Mdouble]: " << std::endl
         << "-hmax [Mdouble]: " << std::endl
             << "Allows to only obtain statistics for specific particles" << std::endl << std::endl
         << "-periodicwalls [bool]: " << std::endl
             << "neglects periodic walls" << std::endl << std::endl
         << "-walls [uint]: " << std::endl
             << "only takes into account the first n walls" << std::endl << std::endl
         << "-verbosity [uint]: " << std::endl
             << "amount of screen output (0 minimal, 1 normal, 2 maximal)" << std::endl << std::endl
         << "-verbose: " << std::endl
             << "identical to -verbosity 2" << std::endl << std::endl
         << "-format [uint]: " << std::endl
             << "???" << std::endl << std::endl
         << "-tmin,tmax,tint [Mdouble]: " << std::endl
             << "redefines the time interval [tmin,tmax] for which statistics are evaluated. Default values from restart file; tint sets tmin to tmax-tint." << std::endl << std::endl
         << "-stepsize [uint]: " << std::endl
             << "to skip some timesteps. Default: 1 " << std::endl << std::endl
         << "-statfilename [std::string]: " << std::endl
             << "to change the name of the output file" << std::endl << std::endl
         << "-timeaverage [bool]: " << std::endl
             << "To average in time (1) or to print statistics for all time steps (0). Default: 1" << std::endl << std::endl
         << "-timevariance [bool]: " << std::endl
             << "to print the time variance; only for time averaged data. Default: ?" << std::endl << std::endl
         << "-gradient: " << std::endl
             << "to plot the first derivative of each statistical value" << std::endl << std::endl
         << "-ignoreFixedParticles: " << std::endl
         << "-StressTypeForFixedParticles: " << std::endl
         << "-mirrorAtDomainBoundary [Mdouble]" << std::endl << std::endl
             << "" << std::endl << std::endl
             << "OUTPUT:" << std::endl << std::endl
             << "1-3: Coordinates " << std::endl
             << "4: Nu " << std::endl
             << "5: Density " << std::endl
             << "6-8: Momentum " << std::endl
             << "9-11: DisplacementMomentum " << std::endl
             << "12-17: Displacement" << std::endl
             << "18-23: MomentumFlux" << std::endl
             << "24-29: DisplacementMomentumFlux" << std::endl
             << "30-32: EnergyFlux " << std::endl
             << "33-41: NormalStress " << std::endl
             << "42-50: TangentialStress " << std::endl
             << "51-53: NormalTraction " << std::endl
             << "54-56: TangentialTraction " << std::endl
             << "57-62: Fabric " << std::endl
             << "63-65: CollisionalHeatFlux " << std::endl
             << "66: Dissipation " << std::endl
             << "67: Potential " << std::endl             
             << "68-70: LocalAngularMomentum " << std::endl
             << "71-79: LocalAngularMomentumFlux " << std::endl
             << "80-88: ContactCoupleStress " << std::endl;
            exit(0);
}

template <StatType T>
void StatisticsVector<T>::set_CG_type(const char* new_) {
    if (!strcmp(new_,"Heaviside")) {
        int dim = 3;
        int numcoeff = 1;
        double heavisidecoeff[] = {1};
        set_Polynomial(heavisidecoeff, numcoeff, dim);
    } else if (!strcmp(new_,"Linear")) {
        int dim = 3;
        int numcoeff = 2;
        double lincoeff[] = {-1,1};
        set_Polynomial(lincoeff, numcoeff, dim);
    } else if (!strcmp(new_,"Lucy")) {
        int dim = 3;
        int numcoeff = 5;
        double lucycoeff[] = {-3,8,-6,0,1};
        set_Polynomial(lucycoeff, numcoeff, dim);
    } else if (!strcmp(new_,"Gaussian")) {
        set_CG_type(Gaussian); return;
    } else if (!strcmp(new_,"HeavisideSphere")) {
        set_CG_type(HeavisideSphere); return;
    } else {
        std::cerr << "Error: unknown CG_type: " << new_ << std::endl;
        exit(-1);
    }
    set_PolynomialName(new_);
    set_CG_type(Polynomial);
}

template <StatType T>
void StatisticsVector<T>::set_CG_type(CG new_) {
    if (new_==Polynomial&&CGPolynomial.get_Order()<0) {std::cerr << "Polynomial is not specified" << std::endl; exit(-1);} 
    CG_type = new_; 
}

template <StatType T>
std::string StatisticsVector<T>::print() {
    std::stringstream ss;
    ss  << "Statistical parameters: name=" << get_name() 
        << ",\n nx=" << nx
        << ", ny=" << ny
        << ", nz=" << nz 
        << ", CG_type=";
    if (get_CG_type()==HeavisideSphere) {
        ss << "HeavisideSphere";
    } else if (get_CG_type()==Gaussian) {
        ss << "Gaussian";
    } else if (get_CG_type()==Polynomial) {
        ss << get_PolynomialName();
    } else { std::cerr << "error in CG_function" << std::endl; exit(-1); }

    ss  << ", w=" << sqrt(w2)
        << ", cutoff=" << cutoff 
        << ",\n tStat=[" << get_tminStat() << "," << get_tmaxStat() << "]"
        << ", xStat=[" << get_xminStat() << "," << get_xmaxStat() << "]"
        << ", yStat=[" << get_yminStat() << "," << get_ymaxStat() << "]"
        << ", zStat=[" << get_zminStat() << "," << get_zmaxStat() << "]"
        << ",\n ignoreFixedParticles=" << ignoreFixedParticles
        << ", StressTypeForFixedParticles=" << StressTypeForFixedParticles
        << ",\n mirrorAtDomainBoundary=" << get_mirrorAtDomainBoundary()
        << ",\n doTimeAverage=" << get_doTimeAverage()
        << ", doVariance=" << get_doVariance()
        << ", doGradient=" << get_doGradient()
        << ", verbosity=" << verbosity
        << std::endl;
    return ss.str();
}

template <StatType T>
std::string StatisticsVector<T>::print_CG() {
    std::stringstream ss;
    ss << "w " << sqrt(get_w2()) 
        << " dim " << get_dim()
        << " domainStat " << get_xminStat() << " " << get_xmaxStat() 
        << " " << get_yminStat() << " " << get_ymaxStat()
        << " " << get_zminStat() << " " << get_zmaxStat()
        << " n " << nx << " " << ny << " " << nz
        << " statType " << statType;
    if (CG_type==Polynomial) ss << " CGPolynomial " << CGPolynomial;
    else ss << " CG_type " << CG_type << " cutoff " << cutoff;
    if (doVariance) ss << " doVariance";
    if (doGradient) ss << " doGradient";
    if (doTimeAverage) {
        ss << " doTimeAverage";
        if (nTimeAverageReset!=-1) ss << " nTimeAverageReset " << nTimeAverageReset;
    }
    if (indSpecies!=-1) ss << " indSpecies " << indSpecies;
    if (rmin!=0) ss << " rmin " << rmin;
    if (rmax!=1e20) ss << " rmax " << rmax;
    if (hmax!=1e20) ss << " hmax " << hmax;
    if (!walls) ss  << " noWalls";
    if (!periodicWalls) ss  << " noPeriodicWalls";
    if (!ignoreFixedParticles) ss   << " includeFixedParticles";
    if (StressTypeForFixedParticles!=1) ss  << " StressTypeForFixedParticles " << StressTypeForFixedParticles;
    return ss.str();
}

template <StatType T>
StatisticsPoint<T> StatisticsVector<T>::average(std::vector<StatisticsPoint<T> > &P) {
    StatisticsPoint<T> avg;
    avg.set_zero();
    for (unsigned int i=0; i<P.size(); i++) avg += P[i];
    avg /= P.size();
    return avg;
}

template <StatType T>
void StatisticsVector<T>::initialize_statistics()
{
    reset_statistics();

    StatisticsVector<T>::set_w2(StatisticsVector<T>::get_w2());
        
    StatisticsVector<T>::set_Positions();
    
    //set tmin and tmax if tint is set
    if (get_tintStat()) {
        set_tminStat(get_t());
        set_tmaxStat(get_tminStat() + get_tintStat());
        std::cout << "tint set, time goes from: " << get_tminStat() << " to: " << get_tmaxStat() << std::endl;
    }

    if (!strcmp(get_stat_filename().c_str(),"")) set_stat_filename(); //set ouput filename to default unless otherwise specified
    if(!open_stat_file(std::fstream::out)) {std::cerr << "Problem opening stat file aborting" <<std::endl; exit(-1);};
        
    get_stat_file() << Points.begin()->write_variable_names() << std::endl;
    get_stat_file() << print_CG() << std::endl;
    
    //std::couts variables
    if (verbosity>1) std::cout << std::endl;
    if (verbosity>0) std::cout << std::endl << print() << std::endl;

}

template <StatType T>
bool StatisticsVector<T>::read_next_from_data_file (unsigned int format) {
    static unsigned int count = 0; 
    
    if (count) for (std::vector<BaseParticle*>::iterator it = getParticleHandler().begin(); it!=getParticleHandler().end(); ++it) {
        (*it)->set_PreviousPosition((*it)->get_Position());
    }

    bool ret_val = MD::read_next_from_data_file (format);

    if (!count) for (std::vector<BaseParticle*>::iterator it = getParticleHandler().begin(); it!=getParticleHandler().end(); ++it) {
        (*it)->set_PreviousPosition((*it)->get_Position());
    }

    count++;
    return ret_val;
}

template <StatType T>
void StatisticsVector<T>::write_statistics()
{ 
    std::cout << "writing statistics for t=" << get_t() << std::endl;

    // write to .stat file
    get_stat_file() << std::left << std::setprecision(8) << std::setw(6)  << get_t()  << std::endl;
    for (unsigned int i=0; i<Points.size(); i++)
        get_stat_file() << Points[i];
    if (get_doGradient()) 
    {
        get_stat_file() << std::left << std::setprecision(8) << std::setw(6)  << get_t()  << std::endl;
        for (unsigned int i=0; i<dx.size(); i++)
            get_stat_file() << dx[i];
        get_stat_file() << std::left << std::setprecision(8) << std::setw(6)  << get_t()  << std::endl;
        for (unsigned int i=0; i<dy.size(); i++)
            get_stat_file() << dy[i];
        get_stat_file() << std::left << std::setprecision(8) << std::setw(6)  << get_t()  << std::endl;
        for (unsigned int i=0; i<dz.size(); i++)
            get_stat_file() << dz[i];
    }
    
}

template <StatType T>
void StatisticsVector<T>::finish_statistics()
{
    if (get_doTimeAverage()) write_time_average_statistics();
    stat_file.close();
}

template <StatType T>
void StatisticsVector<T>::process_statistics(bool usethese)
{
    if (check_current_time_for_statistics()&&usethese)
    {
        if (get_mirrorAtDomainBoundary()) {
            for (unsigned int i=0; i<Points.size(); i++) {
                if (Points[i].mirrorParticle>=0) {
                    //add values to the mirrorParticle
                    Points[Points[i].mirrorParticle]+=Points[i];
                    Points[i].set_zero();
                    if (get_doGradient()) {
                        dx[Points[i].mirrorParticle]+=dx[i];
                        dy[Points[i].mirrorParticle]+=dy[i];
                        dz[Points[i].mirrorParticle]+=dz[i];
                    }
                }
            }
        }
        //~ for (unsigned int i=0; i<Points.size(); i++) 
            //~ Points[i].Displacement /= sqr(Points[i].Density);
        if (get_doTimeAverage()) 
        {
            for (unsigned int i=0; i<timeAverage.size(); i++)
            {
                timeAverage[i] += Points[i];
                if (get_doVariance()) 
                    timeVariance[i] += Points[i].getSquared();
                if (get_doGradient())
                {
                    dxTimeAverage[i] += dx[i];
                    dyTimeAverage[i] += dy[i];
                    dzTimeAverage[i] += dz[i];
                }
            }
            nTimeAverage++;
            if (nTimeAverage==nTimeAverageReset) {
                std::cout << "write time average" << std::endl;
                write_time_average_statistics(); 
                nTimeAverage=0;
                for (unsigned int i=0; i<timeAverage.size(); i++)
                    timeAverage[i].set_zero();
            }
        }
        else 
        {
            write_statistics();
        }
    }
    reset_statistics();
}

template <StatType T>
void StatisticsVector<T>::write_time_average_statistics()
{
        static bool first = true;
        
        if (verbosity) std::cout << std::endl << "averaging " << nTimeAverage << " timesteps " << std::endl;
        if(nTimeAverage==0)
        {
            fprintf(stderr, "\n\n\tERROR :: Can not do TimeAverage of 0 timesteps\n\n");
            //exit(-1);
            return;
        }
        if (first) {
            std::cout << "first" << nTimeAverage << std::endl;
            for (unsigned int i=0; i<timeAverage.size(); i++) { 
                timeAverage[i].firstTimeAverage(nTimeAverage);
                if (get_doVariance()) {
                    timeVariance[i].firstTimeAverage(nTimeAverage);
                    timeVariance[i] -= timeAverage[i].getSquared();
                }
                if (get_doGradient()) {
                    dxTimeAverage[i].firstTimeAverage(nTimeAverage);
                    dyTimeAverage[i].firstTimeAverage(nTimeAverage);
                    dzTimeAverage[i].firstTimeAverage(nTimeAverage);
                }
            }
            first=false;
        } else {
            std::cout << "next" << nTimeAverage << std::endl;
            for (unsigned int i=0; i<timeAverage.size(); i++) { 
                timeAverage[i] /= nTimeAverage;
                if (get_doVariance()) {
                    timeVariance[i] /= nTimeAverage;
                    timeVariance[i] -= timeAverage[i].getSquared();
                }
                if (get_doGradient()) {
                    dxTimeAverage[i] /= nTimeAverage;
                    dyTimeAverage[i] /= nTimeAverage;
                    dzTimeAverage[i] /= nTimeAverage;
                }
            }
        }

        // write average to .stat file
        get_stat_file() << std::left << std::setprecision(8) << std::setw(6) << get_tminStat() << " "  << get_t()  << std::endl;
        for (unsigned int i=0; i<timeAverage.size(); i++)   get_stat_file() << timeAverage[i];
        // write to std::cout
        StatisticsPoint<T> avg = average(timeAverage);
        std::cout << "Averages: " << avg.print() << std::endl;
        
        if (get_doVariance()) {
            // write variance to .stat file
            get_stat_file() << std::left << std::setprecision(8) << std::setw(6) << get_tminStat() << " " << get_t() << " " << std::endl;
            for (unsigned int i=0; i<timeVariance.size(); i++)
                get_stat_file() << timeVariance[i];
            //StatisticsPoint<T> var = average(timeVariance);
        } //end if (get_doVariance)
        
        if (get_doGradient()) {
            // write variance to .stat file
            get_stat_file() << std::left << std::setprecision(8) << std::setw(6) <<  get_tminStat() << " " << get_t() << " " << std::endl;
            for (unsigned int i=0; i<timeAverage.size(); i++) get_stat_file() << dxTimeAverage[i];
            get_stat_file() << std::left << std::setprecision(8) << std::setw(6) <<  get_tminStat() << " " << get_t() << " " << std::endl;
            for (unsigned int i=0; i<timeAverage.size(); i++) get_stat_file() << dyTimeAverage[i];
            get_stat_file() << std::left << std::setprecision(8) << std::setw(6) <<  get_tminStat() << " " << get_t() << " " << std::endl;
            for (unsigned int i=0; i<timeAverage.size(); i++) get_stat_file() << dzTimeAverage[i];
        }
        
        set_tminStat(get_t());

}

template <StatType T>
void StatisticsVector<T>::statistics_from_fstat_and_data()
{
    //This opens the file the data will be recalled from
    set_data_filename();
    open_data_file(std::fstream::in);

    //load first set of timesteps to obtain xdim, ydim, zdim

    //try to read data files; if data.0000 does not exist, try the next, up to .9999
    for (int i=1; i<10000; i++) {
        if(read_next_from_data_file(format)) break;
        set_file_counter(i);
    }

    //set tminstat to smallest time evaluated
    if (get_t()>=get_tminStat()) set_tminStat(get_t()-get_dt());
                
    //This opens the file the force data will be recalled from
    set_fstat_filename();
    open_fstat_file(std::fstream::in);

    //This creates the file statistics will be saved to
    //open_stat_file(std::fstream::out);

    if (verbosity>1) {
        std::cout << "Input from " << get_data_filename() << std::endl;
        std::cout << "Input from " << get_fstat_filename() << std::endl;
        std::cout << "Output to " << get_stat_filename() << std::endl << std::endl;
    }
    
    initialize_statistics();
    
    //std::couts collision time
    Mdouble mass = get_Mass_from_Radius(getLargestParticle()->get_Radius());
    if (verbosity>1) std::cout << "Collision Time " << get_collision_time(mass) 
        << " and restitution coefficient " << get_restitution_coefficient(mass) 
        << " for mass " << mass << std::endl
        << std::endl;
    //std::couts particles
    if (verbosity>2) {MD::print(std::cout,true); std::cout << std::endl;}
    else if (verbosity) {MD::print(std::cout); std::cout << std::endl;}
    
    //this increases the file counter
    if (get_options_data()==2 && get_t()<get_tminStat()) {
        if (verbosity>1) find_next_data_file(get_tminStat());
        else find_next_data_file(get_tminStat(), false);
        for(int i=1;i<get_file_counter();i++) jump_fstat();
        
    }
    
    //Output statistics for each time step 
    std::cout<<"Start statistics"<<std::endl;
    do {
        if (get_t()>get_tmaxStat()) { 
            std::cout << "reached end t="<<std::setprecision (9) <<get_t()<<" tmaxStat=" << std::setprecision (9) <<get_tmaxStat()<< std::endl; break;
        } else if (get_t()>=get_tminStat()) {
            if (verbosity>1) std::cout << "Get statistics for t=" << std::setprecision (9) << get_t() << std::endl;
            else if (verbosity) std::cout << "Get statistics for t=" << std::setprecision (9) << get_t() << " \r";
            if (verbosity>1) std::cout << "#particles="<< getParticleHandler().getNumberOfObjects() << std::endl;
            
            if (!walls) getWallHandler().clear();
            if (!periodicWalls) getBoundaryHandler().clear();
                    
            output_statistics();
            gather_force_statistics_from_fstat_and_data();
            process_statistics(true);

        } else {
            if (verbosity>1) std::cout<<"Jumping statistics t="<<get_t()<<" tminStat()="<<get_tminStat()<<std::endl;
            jump_fstat();
        }
        //if step size>1, skip a number of steps.
        if (get_options_data()!=2) for (unsigned int i=1; i<get_step_size(); i++) {
            //you only get here if you change the step size and you use a single data file
            read_next_from_data_file(format);
            jump_fstat();
        }
    } while (read_next_from_data_file(format));
    //set tmax to largest time evaluated
    if (get_t()<get_tmaxStat()) set_tmaxStat(get_t());

    finish_statistics();
    get_data_file().close();
    get_fstat_file().close();
}

template <StatType T>
void StatisticsVector<T>::auto_setdim () {
    if (getParticleHandler().getNumberOfObjects()<1) return;
    std::vector<BaseParticle*>::iterator it = getParticleHandler().begin();
    set_xmin((*it)->get_Position().X-(*it)->get_Radius());
    set_ymin((*it)->get_Position().Y-(*it)->get_Radius());
    set_zmin((*it)->get_Position().Z-(*it)->get_Radius());
    set_xmax((*it)->get_Position().X+(*it)->get_Radius());
    set_ymax((*it)->get_Position().Y+(*it)->get_Radius());
    set_zmax((*it)->get_Position().Z+(*it)->get_Radius());
    for (std::vector<BaseParticle*>::iterator it = getParticleHandler().begin(); it!=getParticleHandler().end(); ++it) {
        set_xmin(std::min(get_xmin(),(*it)->get_Position().X-(*it)->get_Radius()));
        set_ymin(std::min(get_ymin(),(*it)->get_Position().Y-(*it)->get_Radius()));
        set_zmin(std::min(get_zmin(),(*it)->get_Position().Z-(*it)->get_Radius()));
        set_xmax(std::max(get_xmax(),(*it)->get_Position().X+(*it)->get_Radius()));
        set_ymax(std::max(get_ymax(),(*it)->get_Position().Y+(*it)->get_Radius()));
        set_zmax(std::max(get_zmax(),(*it)->get_Position().Z+(*it)->get_Radius()));
    } //end for all particles
}   

template <StatType T>
bool StatisticsVector<T>::read_next_from_p3p_file () 
{
    //uncomment if you want data file output
    //MD::output_xballs_data();
    
    unsigned int N;
    Mdouble dummy;
    std::string dummyStr;
    
    //read first parameters and check if we reached the end of the file
    
    //Ignoring first line
    getline(p3p_file,dummyStr); 
    std::cout << dummyStr << std::endl; 
    //Reading second line
    p3p_file >> dummy; set_t(dummy);
    p3p_file >> N;
    std::cout << "t= " << get_t() << ": loading " << N << " particles ..." << std::endl;
    getline(p3p_file,dummyStr);
    //std::cout << dummyStr << std::endl;
    //Ignoring third line
    getline(p3p_file,dummyStr);
    std::cout << dummyStr << std::endl;
    if (p3p_file.eof()||p3p_file.peek()==-1) {std::cout << "reached end of p3p file" << std::endl;  return false;}
    
    BaseParticle P0;
    if(getParticleHandler().getNumberOfObjects()<N)
        while(getParticleHandler().getNumberOfObjects()<N)
            getParticleHandler().copyAndAddParticle(P0);
    else
        while(getParticleHandler().getNumberOfObjects()>N)
            getParticleHandler().removeLastParticle();
            
    //Read forth to (N+3)rd line
    Vec3D dummyVec;
    Vec3D position,velocity,angle,angularVelocity;
    int dummyInt;
    for (std::vector<BaseParticle*>::iterator it = getParticleHandler().begin(); it!=getParticleHandler().end(); ++it) {
        p3p_file >> dummyInt; (*it)->set_Index(dummyInt); //ID
        //std::cout << "  ID " << (*it)->get_Index() << std::endl;
        p3p_file >> dummyInt; (*it)->set_species(dummyInt-1); //Group
        p3p_file >> dummy; (*it)->set_Radius(pow(dummy/(4./3.*constants::pi),1./3.)); //Volume
        p3p_file >> dummy; (*it)->set_Mass(dummy); //Mass
        p3p_file 
            >> position
            >> velocity
            >> angularVelocity;
        (*it)->set_Position(position);
        (*it)->set_Velocity(velocity);
        (*it)->set_AngularVelocity(angularVelocity);
        //clean up tangential springs
    } //end for all particles
    //auto_setdim();
    
    getline(p3p_file,dummyStr);
    std::cout << dummyStr << std::endl;
    
    //fix particles, if data_FixedParticles!=0
        
    return true;
}


template <StatType T>
void StatisticsVector<T>::statistics_from_p3()
{
    //This opens the files the data will be recalled from
    std::string p3p_filename = problem_name.str().c_str(); p3p_filename += ".p3p";
    bool opened = open_file (p3p_file, p3p_filename, 1,std::fstream::in);
    std::string p3c_filename = problem_name.str().c_str(); p3c_filename += ".p3c";
    open_file (p3c_file, p3c_filename, 1,std::fstream::in);
    std::string p3w_filename = problem_name.str().c_str(); p3w_filename += ".p3w";
    open_file (p3w_file, p3w_filename, 1,std::fstream::in);
    int version=3;
    //If p3 is not available, open p4 files
    if (!opened) {
        std::cerr << "Warning could not open " << p3p_filename << std::endl; 
        version=4;
        p3p_filename = problem_name.str().c_str(); p3p_filename += ".p4p";
        opened = open_file (p3p_file, p3p_filename, 1,std::fstream::in);
        p3c_filename = problem_name.str().c_str(); p3c_filename += ".p4c";
        opened = open_file (p3c_file, p3c_filename, 1,std::fstream::in);
        p3w_filename = problem_name.str().c_str(); p3w_filename += ".p4w";
        opened = open_file (p3w_file, p3w_filename, 1,std::fstream::in);
        if (!opened) {
            std::cerr << "could not open " << p3p_filename << std::endl; 
            exit(-1);
        }
    }
    std::cout << "opened " << p3p_filename << std::endl;

    //load first set of timesteps to obtain xdim, ydim, zdim

    //try to read data files; if data.0000 does not exist, try the next, up to .9999
    read_next_from_p3p_file();

    //set tminstat to smallest time evaluated
    if (get_t()>=get_tminStat()) set_tminStat(get_t()-get_dt());
    MD::save_restart_data();

    set_data_filename();
    open_data_file(std::fstream::out);
    MD::set_format(14);
    MD::set_dim(3);
    Species[0].set_dim_particle(3);
    MD::output_xballs_data();
    MD::create_xballs_script();
    //exit(-1);
                
    //This creates the file statistics will be saved to
    //open_stat_file(std::fstream::out);

    initialize_statistics();
    
    //std::couts collision time
    Mdouble mass = get_Mass_from_Radius(getLargestParticle()->get_Radius());
    if (verbosity>1) std::cout << "Collision Time " << get_collision_time(mass) 
        << " and restitution coefficient " << get_restitution_coefficient(mass) 
        << " for mass " << mass << std::endl
        << std::endl;
    //std::couts particles
    if (verbosity>2) {MD::print(std::cout,true); std::cout << std::endl;}
    else if (verbosity) {MD::print(std::cout); std::cout << std::endl;}
        
    //Output statistics for each time step 
    std::cout<<"Start statistics"<<std::endl;
    do {
        if (get_t()>get_tmaxStat()) { 
            std::cout << "reached end t="<<std::setprecision (9) <<get_t() <<" tmaxStat=" << std::setprecision (9) <<get_tmaxStat()<< std::endl; 
            break;
        } else if (get_t()>=get_tminStat()) {
            if (verbosity>1) std::cout << "Get statistics for t=" << std::setprecision (9) << get_t() << std::endl;
            else if (verbosity) std::cout << "Get statistics for t=" << std::setprecision (9) << get_t() << " \r";
            if (verbosity>1) std::cout << "#particles="<< getParticleHandler().getNumberOfObjects() << std::endl;
            
            if (!walls) getWallHandler().clear();
            if (!periodicWalls) getBoundaryHandler().clear();
                    
            output_statistics();
            gather_force_statistics_from_p3c(version);
            process_statistics(true);

        } else {
            if (verbosity>1) std::cout 
            << "Jumping statistics t=" << get_t() 
            << " tminStat()=" << get_tminStat() << std::endl;
            jump_p3c();
        }
        //if step size>1, skip a number of steps.
        for (unsigned int i=1; i<get_step_size(); i++) {
            //you only get here if you change the step size and you use a single data file
            read_next_from_p3p_file();
            jump_p3c();
        }
    } while (read_next_from_p3p_file());
    //set tmax to largest time evaluated
    if (get_t()<get_tmaxStat()) set_tmaxStat(get_t());

    finish_statistics();
    p3p_file.close();
    p3c_file.close();
    p3w_file.close();
}

template <StatType T>
void StatisticsVector<T>::jump_p3c()
{
    std::string dummyStr;
    Mdouble N, dummy;
    
    //Ignoring first line
    getline(p3c_file, dummyStr); 
    std::cout << dummyStr << std::endl; 
    //Reading second line
    p3c_file >> dummy; set_t(dummy);
    p3c_file >> N;
    std::cout << "ignoring t= " << get_t() << ": loading " << N << " contacts ..." << std::endl;
    getline(p3c_file,dummyStr);
    //Ignoring third line
    getline(p3c_file,dummyStr);
    std::cout << dummyStr << std::endl;
    //go through N lines in the p3c file
    for (int i=0; i<N; i++) getline(p3c_file,dummyStr);
    
    //Ignoring first line
    getline(p3w_file, dummyStr); 
    std::cout << dummyStr << std::endl; 
    //Reading second line
    p3w_file >> dummy; set_t(dummy);
    p3w_file >> N;
    std::cout << "ignoring t= " << get_t() << ": loading " << N << " wall contacts ..." << std::endl;
    getline(p3w_file,dummyStr);
    //Ignoring third line
    getline(p3w_file,dummyStr);
    //go through N lines in the p3w file
    for (int i=0; i<N; i++) getline(p3w_file,dummyStr);

}

template <StatType T>
void StatisticsVector<T>::gather_force_statistics_from_p3c(int version)
{
    std::string dummyStr;
    Mdouble N, dummy;
    
    //Ignoring first line
    getline(p3c_file, dummyStr); 
    std::cout << dummyStr << std::endl; 
    //Reading second line
    p3c_file >> dummy; set_t(dummy);
    p3c_file >> N;
    std::cout << "t= " << get_t() << ": loading " << N << " contacts ..." << std::endl;
    getline(p3c_file,dummyStr);
    //Ignoring third line
    getline(p3c_file,dummyStr);
    std::cout << dummyStr << std::endl;
    //Checking for end of file
    //if (p3c_file.eof()||p3c_file.peek()==-1) return false;
    
    //set up index list
    int Nmax = 0;
    for (std::vector<BaseParticle*>::iterator it = getParticleHandler().begin(); it!=getParticleHandler().end(); ++it) 
        Nmax = std::max(Nmax,(*it)->get_Index());
    static std::vector<int> index;
    index.resize(Nmax);
    int i=0;
    for (std::vector<BaseParticle*>::iterator it = getParticleHandler().begin(); it!=getParticleHandler().end(); ++it) {
        index[(*it)->get_Index()]=i;
        i++;
    }
    
    // ==> Particledata_4.2to4.25s.csv.p4c <==
    // TIMESTEP CONTACTS
    // 4.2001  31367
    //  P1 P2 CX CY CZ FX FY FZ
    //  9711 9424 0.0549396 0.009705380.149582 4.37283e-05 2.02347e-05 -0.000250662

    // ==> /Users/weinhartt/CarlosLabra/ConfinedCompression_P3P_Spheres.csv.p3c <==
    // TIMESTEP CONTACTS
    // 6.500010000000e-01  21886
    //  P1 P2 FX FY FZ
    //  1126 218 -1.899612000000e-02 6.716474700000e-02 -1.527340000000e-02

    int index1_p3, index2_p3;
    int index1, index2;
    Vec3D Contact, Force;
    Mdouble delta, fdotn, fdott;
    Vec3D P1_P2_normal, P1_P2_tangential;
    Mdouble dist;
    BaseParticle* P1;
    BaseParticle* P2;

    Vec3D contact;
    int counter = 0;
    //go through N lines in the p4c file
    for (int i=0; i<N; i++) {
        counter++;
        if (version==3) {
            p3c_file 
                >> index1_p3
                >> index2_p3
                >> Force;
        } else {
            p3c_file 
                >> index1_p3
                >> index2_p3
                >> Contact
                >> Force;
        }
        index1=index[index1_p3];
        index2=index[index2_p3];
        P1 = getParticleHandler().getObject(index1);
        P2 = getParticleHandler().getObject(index2);
        P1_P2_normal = P1->get_Position() - P2->get_Position();
        dist = GetLength(P1_P2_normal);
        P1_P2_normal /= dist;
        delta = P1->get_Radius() + P2->get_Radius() - dist;
        if (version==3) {
            Contact = P2->get_Position() + P1_P2_normal * (P2->get_Radius()-delta/2);
        }
        fdotn = -Dot(Force,P1_P2_normal)/GetLength2(P1_P2_normal);
        P1_P2_tangential = Force + fdotn * P1_P2_normal; //get tangential force
        fdott = GetLength(P1_P2_tangential);
        if (fdott==0) {
            P1_P2_tangential = Vec3D(0,0,0);
        } else  {   
            P1_P2_tangential /= fdott;
        }
        //Finished reading file
        gather_statistics_collision(index1,index2,Contact,delta,0,fdotn,fdott, P1_P2_normal,-P1_P2_tangential);
        gather_statistics_collision(index2,index1,Contact,delta,0,fdotn,fdott, -P1_P2_normal, P1_P2_tangential);
    }
    
    
    if (verbosity>1) std::cout << "#forces="<< counter << std::endl;
    //Ignoring last line
    getline(p3c_file,dummyStr);
    std::cout << dummyStr << std::endl;
    
    gather_force_statistics_from_p3w(version,index);
}

template <StatType T>
void StatisticsVector<T>::gather_force_statistics_from_p3w(int version, std::vector<int>& index)
{
    std::string dummyStr;
    Mdouble N, dummy;
    
    //Ignoring first line
    getline(p3w_file, dummyStr); 
    std::cout << dummyStr << std::endl; 
    //Reading second line
    p3w_file >> dummy; set_t(dummy);
    p3w_file >> N;
    std::cout << "t= " << get_t() << ": loading " << N << " wall contacts ..." << std::endl;
    getline(p3w_file,dummyStr);
    //Ignoring third line
    getline(p3w_file,dummyStr);
    std::cout << dummyStr << std::endl;
    
    // ==> /Users/weinhartt/CarlosLabra/ConfinedCompression_P3P_Spheres.csv.p3w <==
    // TIMESTEP CONTACTS
    // 6.500010000000e-01  3778
    //  P1 FX FY FZ NX NY NZ (N=Branch)
    //  671 0.000000000000e+00 0.000000000000e+00 1.965970000000e-01 0.000000000000e+00 0.000000000000e+00 -1.991510000000e-03
    //  ==> Particledata_4.2to4.25s.csv.p4w <==
    // TIMESTEP CONTACTS
    // 4.2001  5770
    //  P1 CX CY CZ FX FY FZ
    //   10590 0.00944221 0.0120.130307 -3.14611e-06 -0.000512936 0.000102539
    //   10146 0.0615829 00.162907 2.02968e-05 0.000948522 0.000188616

    int index1_p3;
    int index1, index2;
    Vec3D Contact, Force, Branch;
    Mdouble delta, fdotn, fdott;
    Vec3D P1_P2_normal, P1_P2_tangential;
    Mdouble dist;
    BaseParticle* P1;

    Vec3D contact;
    int counter = 0;
    //go through each line in the fstat file; break when eof or '#' or newline
    for (int i=0; i<N; i++) {
        counter++;
        if (version==3) {
            p3w_file 
                >> index1_p3
                >> Force
                >> Branch;
        } else {
            p3w_file 
                >> index1_p3
                >> Contact
                >> Force;
        }
        index1=index[index1_p3];
        index2=-1;
        P1 = getParticleHandler().getObject(index1);
        if (version==3) {
            //Branch *= -1;
            Contact = P1->get_Position() + Branch;
        } else {
            Branch = Contact - P1->get_Position();
        }
        dist = GetLength(Branch);
        P1_P2_normal = Branch / dist;
        delta = P1->get_Radius() - dist; //check
        fdotn = -Dot(Force,P1_P2_normal)/GetLength2(P1_P2_normal);
        P1_P2_tangential = Force + fdotn * P1_P2_normal; //get tangential force
        fdott = GetLength(P1_P2_tangential);
        if (fdott==0) {
            P1_P2_tangential = Vec3D(0,0,0);
        } else  {   
            P1_P2_tangential /= fdott;
        }
        //Finished reading fstat file
        gather_statistics_collision(index1,index2,Contact,delta,0,fdotn,fdott, P1_P2_normal, -P1_P2_tangential);
    }
    
    if (verbosity>1) std::cout << "#wall forces="<< counter << std::endl;
    //Ignoring last line
    getline(p3w_file,dummyStr);
    std::cout << dummyStr << std::endl;
}





template <StatType T>
void StatisticsVector<T>::set_Positions()
{
    //Dimensions of the domain / n 
    Vec3D diff = Vec3D(
        (get_xmaxStat()-get_xminStat())/nx, 
        (get_ymaxStat()-get_yminStat())/ny,
        (get_zmaxStat()-get_zminStat())/nz);
    //Center of the domain
    Vec3D avg = 0.5 * Vec3D(
        get_xmaxStat()+get_xminStat(), 
        get_ymaxStat()+get_yminStat(), 
        get_zmaxStat()+get_zminStat());
    //Left endpoint of the domain
    Vec3D Min = Vec3D(
        get_xminStat(), 
        get_yminStat(), 
        get_zminStat());

    int num[3]={0,0,0};
    if (get_mirrorAtDomainBoundary()) {
        if (statType==X||statType==XY||statType==XZ||statType==XYZ) {
            //add so many points to domain on both ends
            //todo: the 8 is arbitrary
            num[0] = floor(std::max(get_mirrorAtDomainBoundary(),get_cutoff())/diff.X);
            nx+=2*num[0];
            Min.X-=num[0]*diff.X;
        }
        if (statType==Y||statType==XY||statType==YZ||statType==XYZ) {
            num[1] = floor(std::max(get_mirrorAtDomainBoundary(),get_cutoff())/diff.Y);
            ny+=2*num[1];
            Min.Y-=num[1]*diff.Y;
        }
        if (statType==Z||statType==XZ||statType==XZ||statType==XYZ) {
            num[2] = floor(std::max(get_mirrorAtDomainBoundary(),get_cutoff())/diff.Z);
            nz+=2*num[2];
            Min.Z-=num[2]*diff.Z;
        }
    }

    //Set position of statistics points and set variables to 0
    int N =std::max(nx,1)*std::max(ny,1)*std::max(nz,1);
    Points.resize(N);
    int n = 0;
    for (int i=0; i<std::max(nx,1); i++)
    for (int j=0; j<std::max(ny,1); j++)
    for (int k=0; k<std::max(nz,1); k++) {
        Points[n].set_Position( Vec3D( (nx>1)?(Min.X+diff.X*(i+0.5)):avg.X, 
            (ny>1)?(Min.Y+diff.Y*(j+0.5)):avg.Y, 
            (nz>1)?(Min.Z+diff.Z*(k+0.5)):avg.Z ));
        if (doublePoints) Points[n].set_Position( Points[n].get_Position()
            -Vec3D((i%2)?(0.99*diff.X):0,(j%2)?(0.99*diff.Y):0,(k%2)?(0.99*diff.Z):0) );
        Points[n].set_CG_invvolume();
        Points[n].set_zero();
        if (get_mirrorAtDomainBoundary())
        if ((i<num[0])||(i>nx-num[0])||(j<num[1])||(j>ny-num[1])||(k<num[2])||(k>nz-num[2])) { 
            Points[n].mirrorParticle = n + 
                + ((i<num[0])?(2*(num[0]-i)-1)*ny*nz:0)+((i>nx-num[0])?(2*(nx-num[0]-i)+1)*ny*nz:0)
                + ((j<num[1])?(2*(num[1]-j)-1)*nz:0   )+((j>ny-num[1])?(2*(ny-num[1]-j)+1)*nz:0)
                + ((k<num[2])?(2*(num[2]-k)-1):0      )+((k>nz-num[2])?(2*(nz-num[2]-k)+1):0); 
        }
        n++;
    }
    if (statType==RAZ || statType==RZ  || statType==RA  || statType==AZ || statType==R || statType==A) {
        for (unsigned int k=0; k<Points.size(); k++) {
            Points[k].set_Position(Points[k].get_Position().GetFromCylindricalCoordinates());
        }
    }
    if (get_doGradient()) {
        dx.resize(N);
        dy.resize(N);
        dz.resize(N);
        for (int n=0; n<N; n++) {
            dx[n].set_Position(Points[n].get_Position());
            dx[n].mirrorParticle=Points[n].mirrorParticle;
            dx[n].set_zero();
            dy[n].set_Position(Points[n].get_Position());
            dy[n].mirrorParticle=Points[n].mirrorParticle;
            dy[n].set_zero();
            dz[n].set_Position(Points[n].get_Position());
            dz[n].mirrorParticle=Points[n].mirrorParticle;
            dz[n].set_zero();
        }
    }

    //Set TimeAverage the same way
    if (get_doTimeAverage()) {
        timeAverage.resize(N);
        int n = 0;
        for (int i=0; i<std::max(nx,1); i++)
            for (int j=0; j<std::max(ny,1); j++)
                for (int k=0; k<std::max(nz,1); k++) {
                    timeAverage[n].set_Position(Points[n].get_Position());
                    timeAverage[n].set_zero();
                    timeAverage[n].mirrorParticle=Points[n].mirrorParticle;
                    n++;
                }

        if (get_doVariance()) {
            timeVariance.resize(N);
            for (int n=0; n<N; n++) {
                timeVariance[n].set_Position(timeAverage[n].get_Position());
                timeVariance[n].mirrorParticle=Points[n].mirrorParticle;
                timeVariance[n].set_zero();
            }
        }

        if (get_doGradient()) {
            dxTimeAverage.resize(N);
            dyTimeAverage.resize(N);
            dzTimeAverage.resize(N);
            for (int n=0; n<N; n++) {
                dxTimeAverage[n].set_Position(timeAverage[n].get_Position());
                dxTimeAverage[n].mirrorParticle=Points[n].mirrorParticle;
                dxTimeAverage[n].set_zero();
                dyTimeAverage[n].set_Position(timeAverage[n].get_Position());
                dyTimeAverage[n].mirrorParticle=Points[n].mirrorParticle;
                dyTimeAverage[n].set_zero();
                dzTimeAverage[n].set_Position(timeAverage[n].get_Position());
                dzTimeAverage[n].mirrorParticle=Points[n].mirrorParticle;
                dzTimeAverage[n].set_zero();
            }
        }

    } //end if (get_doTimeAverage()) 
}

template <StatType T>
void StatisticsVector<T>::jump_fstat()
{
    if (get_options_fstat()==2) {
        increase_counter_fstat(std::fstream::in);
        if (verbosity>1) std::cout << "Using fstat file counter: "<<get_file_counter() << std::endl;
    } else {
        //read in first three lines (should start with '#')
        std::string dummy;
        getline (get_fstat_file(), dummy);
        getline (get_fstat_file(), dummy);
        getline (get_fstat_file(), dummy);
        //read in all lines belonging to this timestep
        while ( (get_fstat_file().peek() != -1) && (get_fstat_file().peek() != '#') )
            getline (get_fstat_file(), dummy);
    }
}

template <StatType T>
bool StatisticsVector<T>::satisfiesInclusionCriteria(BaseParticle* P) {
    return P->get_Radius()>rmin 
    && P->get_Radius()<rmax 
    && P->get_Position().Z<hmax;
}

template <StatType T>
void StatisticsVector<T>::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)
{
    Vec3D P1_P2_VelocityAverage, P1_P2_VelocityDifference, P1_P2_Force;
    Mdouble norm_dist;
    if(check_current_time_for_statistics())
    {
        //make sure that for a particle-fixed particle collision the fixed particle is on index2 (switch particles if needed)
        if (getParticleHandler().getObject(index1)->isFixed()) {
            int tmp = index1; index1 = index2; index2 = tmp;
            P1_P2_normal_*=-1;
        }

        if (satisfiesInclusionCriteria(getParticleHandler().getObject(index1)))
        {
            P1_P2_normal=P1_P2_normal_;
            
            if (get_dim()==2) Contact.Z = 0.0;
                
            if (index2<0)
            { //wall-particle collision
                P1_P2_distance = getParticleHandler().getObject(index1)->get_Radius()-delta;
                P2 = Contact;
                P1 = Contact-P1_P2_normal*P1_P2_distance; //note, here we use a minus
                norm_dist=P1_P2_distance;
                P1_P2_VelocityAverage = (getParticleHandler().getObject(index1)->get_Velocity())/2;
                P1_P2_VelocityDifference = (getParticleHandler().getObject(index1)->get_Velocity());
            
                if (StressTypeForFixedParticles==3) {
                    P1_P2_distance += 300;
                    P2 += 300*P1_P2_normal;
                    //this is because wall-particle collisions appear only once in fstat
                    norm_dist=P1_P2_distance;
                }           
            }
            else if (getParticleHandler().getObject(index2)->isFixed()||!satisfiesInclusionCriteria(getParticleHandler().getObject(index2)))
            { //particle-fixed particle collision (external force acts at contact point)
                if (getParticleHandler().getObject(index2)->isFixed()&&StressTypeForFixedParticles==3) { 
                    //infinite extension of stress
                    P1_P2_distance = getParticleHandler().getObject(index1)->get_Radius()+getParticleHandler().getObject(index2)->get_Radius()-delta;
                    P1 = Contact+P1_P2_normal*(0.5*P1_P2_distance); //flowing particle
                    P2 = Contact-P1_P2_normal*(0.5*P1_P2_distance); //fixed particle
                    if (P1_P2_normal.Z<0) { //in this case revert flowing and fixed particle
                        //std::cerr << "Warning: Force normal upwards on base z=" << P1.Z << "." << P2.Z << std::endl;
                        //turning force the other way to achieve sigma=0 at z=inf
                        P1_P2_normal*=-1.0;
                        fdotn *=-1.0;
                    }
                    P1_P2_distance=setInfinitelyLongDistance();
                    P2 = P1 - P1_P2_normal*P1_P2_distance; //move endpoint of force line downwards beyond fixed particles out of the domain
                    if (P1_P2_normal.Z<0) std::cout << P1_P2_distance << " z" << P1.Z << " " << P2.Z  << " f" << fdotn << std::endl;
                } else if (StressTypeForFixedParticles==1|| (!getParticleHandler().getObject(index2)->isFixed()&&StressTypeForFixedParticles==3)) {
                    //add force from flowing particle to contact point to stress
                    P1_P2_distance = getParticleHandler().getObject(index1)->get_Radius()+getParticleHandler().getObject(index2)->get_Radius()-delta;
                    P1 = Contact+P1_P2_normal*(0.5*P1_P2_distance); //1st particle
                    //P2 = Contact; //Contact point;
                    //P1_P2_distance /= 2.0;
                    P1_P2_distance = getParticleHandler().getObject(index1)->get_Radius()-delta/2; //Contact point;
                    P2 = P1 - P1_P2_normal*P1_P2_distance; //Contact point (corrected for polydispersed particles);
                } else {
                    //add force from flowing particle to fixed particle to stress
                    P1_P2_distance = getParticleHandler().getObject(index1)->get_Radius()+getParticleHandler().getObject(index2)->get_Radius()-delta;
                    P1 = Contact+P1_P2_normal*(0.5*P1_P2_distance);
                    P2 = Contact-P1_P2_normal*(0.5*P1_P2_distance);
                }
                norm_dist=P1_P2_distance;//*2.0*Particles[index1].Radius/(Particles[index2].Radius+Particles[index1].Radius);
                //This is the length of the contact in particle 1.
                P1_P2_VelocityAverage = getParticleHandler().getObject(index2)->get_Velocity()/2;
                P1_P2_VelocityDifference = -getParticleHandler().getObject(index2)->get_Velocity();
            }
            else
            {//particle-particle collision

                if (getParticleHandler().getObject(index2)->isFixed()) std::cout << "ERROR" << std::endl;
                P1_P2_distance = getParticleHandler().getObject(index1)->get_Radius()+getParticleHandler().getObject(index2)->get_Radius()-delta;

                P1 = Contact+P1_P2_normal*(0.5*P1_P2_distance);
                P2 = Contact-P1_P2_normal*(0.5*P1_P2_distance);
                //This is the length of the contact in particle 1.
                norm_dist=P1_P2_distance*getParticleHandler().getObject(index1)->get_Radius()/(getParticleHandler().getObject(index2)->get_Radius()+getParticleHandler().getObject(index1)->get_Radius());
                P1_P2_VelocityAverage = (getParticleHandler().getObject(index1)->get_Velocity()+getParticleHandler().getObject(index2)->get_Velocity())/2;
                P1_P2_VelocityDifference = (getParticleHandler().getObject(index1)->get_Velocity()-getParticleHandler().getObject(index2)->get_Velocity());
            }
            //Particles[max(index1,index2)] is chosen so that it is the non-wall, non-fixed particle
            P1_P2_Fabric = SymmetrizedDyadic(P1_P2_normal, P1_P2_normal) * getParticleHandler().getObject(index1)->get_Volume(get_Species());
            //fix
            //~ P1_P2_Fabric = SymmetrizedDyadic(P1_P2_normal, P1_P2_normal) * getParticleHandler().getObject(max(index1,index2))->get_Volume(get_Species());
            P1_P2_NormalStress = Dyadic(P1_P2_normal,P1_P2_normal) * (fdotn*norm_dist);
            P1_P2_TangentialStress = Dyadic(P1_P2_tangential,P1_P2_normal) * (fdott*norm_dist);
            P1_P2_NormalTraction     = P1_P2_normal     * fdotn;
            P1_P2_TangentialTraction = P1_P2_tangential * fdott;
            P1_P2_ContactCoupleStress = Dyadic(P1_P2_NormalTraction+P1_P2_TangentialTraction, -0.5*P1_P2_normal*norm_dist);
            P1_P2_Contact = Contact;
            P1_P2_Potential = (get_k()*sqr(delta)+get_kt()*sqr(ctheta))/2;
            P1_P2_Force = fdotn*P1_P2_normal+fdott*P1_P2_tangential;
            P1_P2_Dissipation = Dot(P1_P2_VelocityDifference,P1_P2_Force)/2;
            P1_P2_CollisionalHeatFlux = P1_P2_normal * (P1_P2_distance/2 * Dot(P1_P2_VelocityAverage,fdotn*P1_P2_normal+fdott*P1_P2_tangential)) + P1_P2_Potential * P1_P2_VelocityAverage;
            //Now P1_P2 distance and normal are in non-averaged directions - This is the projection of the distance on to the CG directions. Required for the later statistics hence done now, not before.
            P1_P2_distance *= sqrt(Points[0].dot(P1_P2_normal,P1_P2_normal));
            
            //If the normal to the contact and the averaging direction are EXACTLY+- parallel then set the normal to zero. Otherwise you would end up with nan, which is not good.
            if (P1_P2_distance>1e-12)
                P1_P2_normal = (P1-P2)/P1_P2_distance;
            else
                P1_P2_normal = Vec3D(0.0,0.0,0.0);
            //P1_P2_normal = (P1-P2)/P1_P2_distance;
            
            //if 2nd particle is wall, fixed particle or not-included particle
            if (index2<0 || getParticleHandler().getObject(index2)->isFixed() || !satisfiesInclusionCriteria(getParticleHandler().getObject(index2)))
            { 
                // << "/" << P2 << std::endl;
                if (StressTypeForFixedParticles!=0) { //only if wall - flow contact adds anything to the stress
                    evaluate_force_statistics();
                }
                if (StressTypeForFixedParticles!=3 || !(index2<0 || getParticleHandler().getObject(index2)->isFixed())) { //only if stress is not infinitely extended
                    Vec3D P;
                    if (StressTypeForFixedParticles==0) {
                        P=P1; //Traction at flow particle
                    } else {
                        P=P2; //Traction at contact point, resp wall particle
                    }           
                    evaluate_wall_force_statistics(P);
                }
            }
            else
            {
                evaluate_force_statistics();
            }
        }
    }
}

template <StatType T>
void StatisticsVector<T>::gather_force_statistics_from_fstat_and_data()
{
    //read in first three lines (should start with '#')
    std::string dummy;
    if (get_options_fstat()==2) increase_counter_fstat(std::fstream::in);
    getline (get_fstat_file(), dummy);
    getline (get_fstat_file(), dummy);
    getline (get_fstat_file(), dummy);
        
            
    static Mdouble time;
    static int index1, index2;
    static Vec3D Contact;
    static Mdouble delta, ctheta, fdotn, fdott;
    static Vec3D P1_P2_normal_, P1_P2_tangential;
    static BaseParticle PI;
    static BaseParticle PJ;

    int counter = 0;
    //go through each line in the fstat file; break when eof or '#' or newline
    while ( (get_fstat_file().peek() != -1) && (get_fstat_file().peek() != '#') ) { //(!get_fstat_file().eof())
        counter++;
        /* # 1: time
         # 2: particle Number i
         # 3: contact partner j (particles >= 0, walls < 0)
         # 4: x-position \
         # 5: y-position  > of the contact point (I hope)
         # 6: z-position /
         # 7: delta = overlap at the contact
         # 8: ctheta = length of the tangential spring
         # 9: P1_P2_normal force |f^n|
         # 10: remaining (tangential) force |f^t|=|f-f^n|
         # 11-13: P1_P2_normal unit vector nx, ny, nz
         # 14-16: tangential unit vector tx, ty, tz
        */
        get_fstat_file() 
            >> time
            >> index1
            >> index2
            >> Contact
            >> delta
            >> ctheta
            >> fdotn
            >> fdott
            >> P1_P2_normal_
            >> P1_P2_tangential;
        get_fstat_file().ignore(256,'\n');
        //Finished reading fstat file
        gather_statistics_collision(index1,index2,Contact,delta,ctheta,fdotn,fdott, P1_P2_normal_, P1_P2_tangential);

    }
    if (verbosity>1) std::cout << "#forces="<< counter << std::endl;
}

template <StatType T>
void StatisticsVector<T>::evaluate_force_statistics(int wp)
{
    if(periodicWalls)
        for (unsigned int k=wp; k<getBoundaryHandler().getNumberOfObjects(); k++)
        {
            PeriodicBoundary* b0=dynamic_cast<PeriodicBoundary*>(getBoundaryHandler().getObject(k));
            if(b0)
            {
                b0->distance(P1);
                b0->shift_positions(P1,P2);
                evaluate_force_statistics(k+1);
                b0->shift_positions(P1,P2);
            }
        }
        
    //evaluate fields that only depend on CParticle parameters
    Mdouble psi; //Course graining integral
    Vec3D rpsi=Vec3D(0,0,0); 
    static unsigned int counter = 0;
    for (unsigned int i=0; i<Points.size(); i++) {
        psi = Points[i].CG_integral(P1, P2, P1_P2_normal, P1_P2_distance, rpsi);
        if (psi) {
            counter++;
            if (!std::isfinite(psi)) {
                std::cerr << "error: psi =" << psi << " is infinite for P1=" << P1 << ", P2=" << P2 << ", t=" << get_t() << std::endl;
                psi=0;
            }
            Points[i].ContactCoupleStress += (Cross(P1_P2_Contact*psi,P1_P2_ContactCoupleStress) - Cross(rpsi,P1_P2_ContactCoupleStress))*(-0.5);
            Points[i].NormalStress += P1_P2_NormalStress * psi;
            Points[i].TangentialStress += P1_P2_TangentialStress * psi;
            Points[i].Fabric += P1_P2_Fabric * psi;
            Points[i].CollisionalHeatFlux += P1_P2_CollisionalHeatFlux * psi;
            Points[i].Dissipation += P1_P2_Dissipation * psi;
            Points[i].Potential += P1_P2_Potential * psi;
            if (get_doGradient()) {
                Vec3D d_psi = Points[i].CG_integral_gradient(P1, P2, P1_P2_normal, P1_P2_distance);

                dx[i].NormalStress += P1_P2_NormalStress * d_psi.X;
                dx[i].TangentialStress += P1_P2_TangentialStress * d_psi.X;
                dx[i].Fabric += P1_P2_Fabric * d_psi.X;
                dx[i].CollisionalHeatFlux += P1_P2_CollisionalHeatFlux * d_psi.X;
                dx[i].Dissipation += P1_P2_Dissipation * d_psi.X;
                dx[i].Potential += P1_P2_Potential * d_psi.X;

                dy[i].NormalStress += P1_P2_NormalStress * d_psi.Y;
                dy[i].TangentialStress += P1_P2_TangentialStress * d_psi.Y;
                dy[i].Fabric += P1_P2_Fabric * d_psi.Y;
                dy[i].CollisionalHeatFlux += P1_P2_CollisionalHeatFlux * d_psi.Y;
                dy[i].Dissipation += P1_P2_Dissipation * d_psi.Y;
                dy[i].Potential += P1_P2_Potential * d_psi.Y;

                dz[i].NormalStress += P1_P2_NormalStress * d_psi.Z;
                dz[i].TangentialStress += P1_P2_TangentialStress * d_psi.Z;
                dz[i].Fabric += P1_P2_Fabric * d_psi.Z;
                dz[i].CollisionalHeatFlux += P1_P2_CollisionalHeatFlux * d_psi.Z;
                dz[i].Dissipation += P1_P2_Dissipation * d_psi.Z;
                dz[i].Potential += P1_P2_Potential * d_psi.Z;
            } //end if get_doGradient
        } //end if phi
    } // end forall Points
    //std::cout << "NCeval" << counter << std::endl;
}
        

template <StatType T>
void StatisticsVector<T>::evaluate_wall_force_statistics(Vec3D P, int wp)
{
    if(periodicWalls)
        for (unsigned int k=wp; k<getBoundaryHandler().getNumberOfObjects(); k++)
        {
            PeriodicBoundary* b0=dynamic_cast<PeriodicBoundary*>(getBoundaryHandler().getObject(k));
            if(b0)
            {
                b0->distance(P1);
                b0->shift_positions(P1,P2);
                evaluate_force_statistics(k+1);
                b0->shift_positions(P1,P2);
            }
        }

    //evaluate fields that only depend on CParticle parameters
    Mdouble phi; //Course graining integral
    for (unsigned int i=0; i<Points.size(); i++) {
        phi = Points[i].CG_function(P);
        if (phi) {
            Points[i].NormalTraction += P1_P2_NormalTraction * phi;
            Points[i].TangentialTraction += P1_P2_TangentialTraction * phi;
            if (get_doGradient()) {
                Vec3D d_phi = Points[i].CG_gradient(P, phi);
                dx[i].NormalTraction += P1_P2_NormalTraction * d_phi.X;
                dy[i].NormalTraction += P1_P2_NormalTraction * d_phi.Y;
                dz[i].NormalTraction += P1_P2_NormalTraction * d_phi.Z;
                dx[i].TangentialTraction += P1_P2_TangentialTraction * d_phi.X;
                dy[i].TangentialTraction += P1_P2_TangentialTraction * d_phi.Y;
                dz[i].TangentialTraction += P1_P2_TangentialTraction * d_phi.Z;
            } //end if get_doGradient
        } //end if phi
    } // end forall Points

}

template <StatType T>
void StatisticsVector<T>::output_statistics() {
    if(check_current_time_for_statistics())
    {
        for (unsigned int i=0; i<Points.size(); i++)
            Points[i].set_CG_invvolume();
            
        for (std::vector<BaseParticle*>::iterator P = getParticleHandler().begin(); P!=getParticleHandler().end(); ++P) {
            //unfix particles (turn off to ignore fixed particles)
            if ((!ignoreFixedParticles) && (*P)->isFixed()) {
                (*P)->unfix(get_Species());
            }
            //ignore fixed particles

            if (satisfiesInclusionCriteria(*P) && !(*P)->isFixed()) {
                evaluate_particle_statistics(P);
            }

        }
    }
}

template <StatType T>
void StatisticsVector<T>::evaluate_particle_statistics(std::vector<BaseParticle*>::iterator P, int wp) {
    bool doDisplacement = false;
    
    if (periodicWalls)
        for (unsigned int k=wp; k<getBoundaryHandler().getNumberOfObjects(); k++)
        {
            PeriodicBoundary* b0=dynamic_cast<PeriodicBoundary*>(getBoundaryHandler().getObject(k));
            if(b0)
            {
                Mdouble dist = b0->distance(P1);
                if (sqr(dist-(*P)->get_Radius())<9.0*get_w2())
                {
                    b0->shift_position(*P);
                    evaluate_particle_statistics(P, k+1);
                    b0->shift_position(*P);
                }               
            }
        }

    //evaluate fields that only depend on CParticle parameters
    Mdouble phi; //Course graining function
    Vec3D V = Vec3D(0,0,0);
    if (VelocityProfile.size()) V = getVelocityProfile ((*P)->get_Position());

    for (unsigned int i=0; i<Points.size(); i++) {
        phi = Points[i].CG_function((*P)->get_Position());
        if (phi) {
            if (!std::isfinite(phi)) {
                std::cout 
                << "t =" << MD::get_t()
                << "error: phi =" << phi 
                << " is infinite for P=" << (*P)->get_Position()
                << ", Point=" << Points[i].get_Position()
                << std::endl;
                phi=0;
            }
            Points[i].Nu += (*P)->get_Volume(get_Species()) * phi;
            Points[i].Density += (*P)->get_Mass() * phi;
            if (VelocityProfile.size()) {
                Points[i].Momentum += ((*P)->get_Velocity()-V) * ((*P)->get_Mass() * phi);
                Points[i].MomentumFlux += SymmetrizedDyadic((*P)->get_Velocity()-V, (*P)->get_Velocity()-V) * ((*P)->get_Mass() * phi);
            } else {
                Points[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * phi);
                Points[i].MomentumFlux += SymmetrizedDyadic((*P)->get_Velocity(), (*P)->get_Velocity()) * ((*P)->get_Mass() * phi);
            }
            Points[i].DisplacementMomentum += (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()) * ((*P)->get_Mass() * phi);
            Points[i].DisplacementMomentumFlux += SymmetrizedDyadic((*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()), (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount())) * ((*P)->get_Mass() * phi);
            Points[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * phi);
            
            Vec3D LocalAngularMomentum = Cross(Points[0].clearAveragedDirections((*P)->get_Position()-Points[i].get_Position()), (*P)->get_Velocity()) * (*P)->get_Mass() 
                + (*P)->get_AngularVelocity() * (*P)->get_Inertia();
            Points[i].LocalAngularMomentum += phi * LocalAngularMomentum;
            Points[i].LocalAngularMomentumFlux += Dyadic(LocalAngularMomentum, (*P)->get_Velocity() * phi);

            //Note: Displacement is only computed directly if gradients are cmputed
            if (get_doGradient()) {
                
                Vec3D d_phi = Points[i].CG_gradient((*P)->get_Position(), phi);
                
                if (doDisplacement) {
                    //begin: Linear displacement, \f$1/(2 \rho^2) \sum_{pq} m_p m_q \phi_q *(v_{pqa} \partial_b \phi_p + v_{pqb} \partial_a \phi_p) \f$ 
                    for (std::vector<BaseParticle*>::iterator Q = getParticleHandler().begin(); Q!=getParticleHandler().end(); ++Q) {
                        double phiQ = Points[i].CG_function((*Q)->get_Position()); //Course graining function
                        Points[i].Displacement += SymmetrizedDyadic((*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()) - (*Q)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()), d_phi) 
                         * ((*P)->get_Mass() * (*Q)->get_Mass() * phiQ);
                    }
                    //end:Displacement
                }
                
                dx[i].Nu += (*P)->get_Volume(get_Species()) * d_phi.X;
                dx[i].Density += (*P)->get_Mass() * d_phi.X;
                dx[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * d_phi.X);
                dx[i].DisplacementMomentum += (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()) * ((*P)->get_Mass() * d_phi.X);
                dx[i].MomentumFlux += SymmetrizedDyadic((*P)->get_Velocity(), (*P)->get_Velocity()) * ((*P)->get_Mass() * d_phi.X);
                dx[i].DisplacementMomentumFlux += SymmetrizedDyadic((*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()), (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount())) * ((*P)->get_Mass() * d_phi.X);
                dx[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * d_phi.X);

                dy[i].Nu += (*P)->get_Volume(get_Species()) * d_phi.Y;
                dy[i].Density += (*P)->get_Mass() * d_phi.Y;
                dy[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * d_phi.Y);
                dy[i].DisplacementMomentum += (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()) * ((*P)->get_Mass() * d_phi.Y);
                dy[i].MomentumFlux += SymmetrizedDyadic((*P)->get_Velocity(), (*P)->get_Velocity()) * ((*P)->get_Mass() * d_phi.Y);
                dy[i].DisplacementMomentumFlux += SymmetrizedDyadic((*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()), (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount())) * ((*P)->get_Mass() * d_phi.Y);
                dy[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * d_phi.Y);

                dz[i].Nu += (*P)->get_Volume(get_Species()) * d_phi.Z;
                dz[i].Density += (*P)->get_Mass() * d_phi.Z;
                dz[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * d_phi.Z);
                dz[i].DisplacementMomentum += (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()) * ((*P)->get_Mass() * d_phi.Z);
                dz[i].MomentumFlux += SymmetrizedDyadic((*P)->get_Velocity(), (*P)->get_Velocity()) * ((*P)->get_Mass() * d_phi.Z);
                dz[i].DisplacementMomentumFlux += SymmetrizedDyadic((*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()), (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount())) * ((*P)->get_Mass() * d_phi.Z);
                dz[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * d_phi.Z);
            } //end if get_doGradient
        } //end if phi
    } // end forall Points
}

//template class StatisticsPoint<XYZ>;

//Templated functions

template<> void StatisticsVector<XY>::set_nz(int new_ UNUSED) { nz = 1; }
template<> void StatisticsVector<XZ>::set_ny(int new_ UNUSED) { ny = 1; }
template<> void StatisticsVector<YZ>::set_nx(int new_ UNUSED) { nx = 1; }
template<> void StatisticsVector<X>::set_ny(int new_ UNUSED) { ny = 1; }
template<> void StatisticsVector<X>::set_nz(int new_ UNUSED) { nz = 1; }
template<> void StatisticsVector<Y>::set_nx(int new_ UNUSED) { nx = 1; }
template<> void StatisticsVector<Y>::set_nz(int new_ UNUSED) { nz = 1; }
template<> void StatisticsVector<Z>::set_nx(int new_ UNUSED) { nx = 1; }
template<> void StatisticsVector<Z>::set_ny(int new_ UNUSED) { ny = 1; }
template<> void StatisticsVector<O>::set_nx(int new_ UNUSED) { nx = 1; }
template<> void StatisticsVector<O>::set_ny(int new_ UNUSED) { ny = 1; }
template<> void StatisticsVector<O>::set_nz(int new_ UNUSED) { nz = 1; }

template<> void StatisticsVector<RAZ>::set_statType() { statType=RAZ; }
template<> void StatisticsVector<RZ>::set_statType() { statType=RZ; }
template<> void StatisticsVector<AZ>::set_statType() { statType=AZ; }
template<> void StatisticsVector<RA>::set_statType() { statType=RA; }
template<> void StatisticsVector<A>::set_statType() { statType=A; }
template<> void StatisticsVector<R>::set_statType() { statType=R; }
template<> void StatisticsVector<XYZ>::set_statType() { statType=XYZ; }
template<> void StatisticsVector<XY>::set_statType() { statType=XY; }
template<> void StatisticsVector<XZ>::set_statType() { statType=XZ; }
template<> void StatisticsVector<YZ>::set_statType() { statType=YZ; }
template<> void StatisticsVector<X>::set_statType() { statType=X; }
template<> void StatisticsVector<Y>::set_statType() { statType=Y; }
template<> void StatisticsVector<Z>::set_statType() { statType=Z; }
template<> void StatisticsVector<O>::set_statType() { statType=O; }

template<> Mdouble StatisticsVector<RAZ>::setInfinitelyLongDistance() {
    return std::min(std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
                   std::max((get_xmaxStat()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X)),
                   std::max((get_ymaxStat()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y));
}
template<> Mdouble StatisticsVector<RA>::setInfinitelyLongDistance() {
    return std::min(std::max((get_ymaxStat()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y),
               std::max((get_xmaxStat()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X));
}
template<> Mdouble StatisticsVector<RZ>::setInfinitelyLongDistance() {
    return std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
               std::max((get_xmaxStat()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X));
}
template<> Mdouble StatisticsVector<AZ>::setInfinitelyLongDistance() {
    return std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
               std::max((get_ymaxStat()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y));
}
template<> Mdouble StatisticsVector<R>::setInfinitelyLongDistance() {
    return std::max((get_xmaxStat()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X);
}
template<> Mdouble StatisticsVector<A>::setInfinitelyLongDistance() {
    return std::max((get_ymaxStat()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y);
}


template<> Mdouble StatisticsVector<XYZ>::setInfinitelyLongDistance() {
    return std::min(std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
                   std::max((get_xmaxStat()-P1.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P1.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X)),
                   std::max((get_ymaxStat()-P1.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P1.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y));
}
template<> Mdouble StatisticsVector<XY>::setInfinitelyLongDistance() {
    return std::min(std::max((get_ymaxStat()-P1.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P1.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y),
               std::max((get_xmaxStat()-P1.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P1.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X));
}
template<> Mdouble StatisticsVector<XZ>::setInfinitelyLongDistance() {
    return std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
               std::max((get_xmaxStat()-P1.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P1.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X));
}
template<> Mdouble StatisticsVector<YZ>::setInfinitelyLongDistance() {
    return std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
               std::max((get_ymaxStat()-P1.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P1.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y));
}
template<> Mdouble StatisticsVector<X>::setInfinitelyLongDistance() {
    return std::max((get_xmaxStat()-P1.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P1.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X);
}
template<> Mdouble StatisticsVector<Y>::setInfinitelyLongDistance() {
    return std::max((get_ymaxStat()-P1.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P1.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y);
}
template<> Mdouble StatisticsVector<Z>::setInfinitelyLongDistance() {
    return fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z);
}
template<> double StatisticsVector<O>::setInfinitelyLongDistance() {
    return 0;
}

• StatisticsVector.hcc
• Generated on Thu Nov 10 2016 11:54:53 for MercuryDPM by 1.8.7