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.


#include "Particles/BaseParticle.h"
#include "Boundaries/PeriodicBoundary.h"
#include "Boundaries/AngledPeriodicBoundary.h"
#include "Species/LinearViscoelasticSpecies.h"

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>::setCGWidth2(Mdouble new_)
{
    // set w2
    if (new_ > 0.0)
        w2 = new_;
    else
        w2 = mathsFunc::square(w_over_rmax * particleHandler.getLargestParticle()->getRadius());
    //(2.0*getXMax()/nx);
    
    // set cutoff radius
    if (getCGShape() == HeavisideSphere || getCGShape() == Polynomial)
    {
        cutoff = sqrt(w2);
    }
    else if (getCGShape() == Gaussian)
    {
        if (getSuperExact())
            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 (getDoGradient())
        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 = particleHandler.begin(); it != particleHandler.end(); ++it)
    {
        (*it)->setPreviousPosition((*it)->getPosition());
    }
}

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;
    setCGTimeAveragingInterval(0);
    setCGTimeMin(-1e20);
    setTimeMaxStat(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;
    setSuperExact(false);
    VelocityProfile.resize(0); //equivalent to do not use velocityprofile
    setStepSize(1);
    
    // time averaging
    setDoTimeAverage(true);
    nTimeAverage = 0;
    nTimeAverageReset = -1;
    //calculate variance
    setDoVariance(false);
    //calculate gradient
    setDoGradient(false);
    doDoublePoints = false;
    
    // additional stuff
    statFile.setFileType(FileType::ONE_FILE);
}

template<StatType T>
StatisticsVector<T>::StatisticsVector(StatisticsVector& other)
        : DPMBase()
{
    //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
    setName(name.c_str());

    if (!strcmp(getName().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;
        if (readRestartFile()==0)
        {
            std::cout << "no restart file given: using default density 1" << std::endl;
            speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());

        }
        setAppend(false);
        //~ setTimeMaxStat(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(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);
    int 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);
    
    setNX(nx);
    setNY(ny);
    setNZ(nz);
    setXMinStat(xmin);
    setYMinStat(ymin);
    setZMinStat(zmin);
    setXMaxStat(xmax);
    setYMaxStat(ymax);
    setZMaxStat(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((Position.X - VelocityProfile_Min.X) / VelocityProfile_D.X, &i);
    m = modf((Position.Y - VelocityProfile_Min.Y) / VelocityProfile_D.Y, &j);
    n = modf((Position.Z - VelocityProfile_Min.Z) / VelocityProfile_D.Z, &k);
    if (nx < 1)
        i = 0.0;
    if (ny < 1)
        j = 0.0;
    if (nz < 1)
        k = 0.0;
    int ix = static_cast<int>((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(int argc, char *argv[])
{
    //check for options (standardly '-option argument')
    for (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"))
            {
                setCGShape(HeavisideSphere);
            }
            else if (!strcmp(argv[i + 1], "Gaussian"))
            {
                setCGShape(Gaussian);
            }
            else
            {
                setCGShape(argv[i + 1]);
            }
        }
        else if (!strcmp(argv[i], "-w"))
        {
            double old = getCGWidth();
            setCGWidth(atof(argv[i + 1]));
            std::cout << " w changed from " << old << " to " << getCGWidth() << 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"))
        {
            setN(atoi(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-nx"))
        {
            setNX(atoi(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-ny"))
        {
            setNY(atoi(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-nz"))
        {
            setNZ(atoi(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-x"))
        {
            setXMinStat(atof(argv[i + 1]));
            setXMaxStat(atof(argv[i + 2]));
            i++; //requires two arguments
        }
        else if (!strcmp(argv[i], "-y"))
        {
            setYMinStat(atof(argv[i + 1]));
            setYMaxStat(atof(argv[i + 2]));
            i++; //requires two arguments
        }
        else if (!strcmp(argv[i], "-z"))
        {
            setZMinStat(atof(argv[i + 1]));
            setZMaxStat(atof(argv[i + 2]));
            i++; //requires two arguments
        }
        else if (!strcmp(argv[i], "-positions"))
        {
            loadPositions(argv[i+1]);
        }
        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"))
        {
            setDoPeriodicWalls(atoi(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-walls"))
        {
            setCGWidthalls(atoi(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-verbosity"))
        {
            setVerbosityLevel(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], "-doDoublePoints"))
        {
            if (argc <= i + 1 || argv[i + 1][0] == '-')
            {
                setDoDoublePoints(true);
                i--; //no argument
            }
            else
                setDoDoublePoints(atoi(argv[i + 1]));
            std::cout << "set doDoublePoints=" << getDoDoublePoints() << std::endl;
        }
        else if (!strcmp(argv[i], "-w_over_rmax"))
        {
            setCGWidth_over_rmax(atof(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-tmin"))
        {
            setCGTimeMin(atof(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-tmax"))
        {
            setTimeMaxStat(atof(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-tint"))
        {
            setCGTimeAveragingInterval(atof(argv[i + 1]));
            //~ } else if (!strcmp(argv[i],"-timesteps")) {
            //~ setCGTimeAveragingInterval(atof(argv[i+1])*getTimeStep());
        }
        else if (!strcmp(argv[i], "-stepsize"))
        {
            int stepSize = atoi(argv[i + 1]);
            setStepSize(stepSize);
        }
        else if (!strcmp(argv[i], "-loaddatafile"))
        {
            readDataFile(argv[i + 1]);
            setCGTimeMin(getTime());
            setTimeMaxStat(getTimeMax());
        }
        else if (!strcmp(argv[i], "-statfilename") | !strcmp(argv[i], "-o"))
        {
            statFile.setName(argv[i + 1]);
            // setName(argv[i+1]);
        }
        else if (!strcmp(argv[i], "-timeaverage"))
        {
            setDoTimeAverage(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] == '-')
            {
                setDoGradient(true);
                i--; //no argument
            }
            else
                setDoGradient(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] == '-')
            {
                setDoVariance(true);
                i--; //no argument
            }
            else
                setDoVariance(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] == '-')
            {
                setSuperExact(true);
                i--; //no argument
            }
            else
                setSuperExact(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] == '-')
            {
                setDoIgnoreFixedParticles(true);
                i--; //no argument
            }
            else
                setDoIgnoreFixedParticles(atoi(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-StressTypeForFixedParticles"))
        {
            StressTypeForFixedParticles = atoi(argv[i + 1]);
        }
        else if (!strcmp(argv[i], "-mirrorAtDomainBoundary"))
        {
            setMirrorAtDomainBoundary(atof(argv[i + 1]));
        }
        else if (!strcmp(argv[i], "-stattype") || !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>::setCGShape(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"))
    {
        setCGShape(Gaussian);
        return;
    }
    else if (!strcmp(new_, "HeavisideSphere"))
    {
        setCGShape(HeavisideSphere);
        return;
    }
    else
    {
        std::cerr << "Error: unknown CG_type: " << new_ << std::endl;
        exit(-1);
    }
    setPolynomialName(new_);
    setCGShape(Polynomial);
}

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

template<StatType T>
std::string StatisticsVector<T>::printStat()
{
    std::stringstream ss;
    ss << "Statistical parameters: name=" << getName()
            << ",\n nx=" << nx
            << ", ny=" << ny
            << ", nz=" << nz
            << ", CG_type=";
    if (getCGShape() == HeavisideSphere)
    {
        ss << "HeavisideSphere";
    }
    else if (getCGShape() == Gaussian)
    {
        ss << "Gaussian";
    }
    else if (getCGShape() == Polynomial)
    {
        ss << getPolynomialName();
    }
    else
    {
        std::cerr << "error in CG_function" << std::endl;
        exit(-1);
    }
    
    ss << ", w=" << sqrt(w2)
            << ", cutoff=" << cutoff
            << ",\n tStat=[" << getCGTimeMin() << "," << getTimeMaxStat() << "]"
            << ", xStat=[" << getXMinStat() << "," << getXMaxStat() << "]"
            << ", yStat=[" << getYMinStat() << "," << getYMaxStat() << "]"
            << ", zStat=[" << getZMinStat() << "," << getZMaxStat() << "]"
            << ",\n ignoreFixedParticles=" << ignoreFixedParticles
            << ", StressTypeForFixedParticles=" << StressTypeForFixedParticles
            << ",\n mirrorAtDomainBoundary=" << getMirrorAtDomainBoundary()
            << ",\n doTimeAverage=" << getDoTimeAverage()
            << ", doVariance=" << getDoVariance()
            << ", doGradient=" << getDoGradient()
            << ", verbosity=" << verbosity
            << std::endl;
    return ss.str();
}

template<StatType T>
std::string StatisticsVector<T>::print_CG()
{
    std::stringstream ss;
    ss << "w " << sqrt(getCGWidthSquared())
            << " dim " << getSystemDimensions()
            << " domainStat " << getXMinStat() << " " << getXMaxStat()
            << " " << getYMinStat() << " " << getYMaxStat()
            << " " << getZMinStat() << " " << getZMaxStat()
            << " 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>::initialiseStatistics()
{
    reset_statistics();
    
    StatisticsVector<T>::setCGWidth2(StatisticsVector<T>::getCGWidthSquared());
    
    StatisticsVector<T>::setPositions();
    
    //set tmin and tmax if tint is set
    if (getCGTimeAveragingInterval() != 0.0)
    {
        setCGTimeMin(getTime());
        setTimeMaxStat(getCGTimeMin() + getCGTimeAveragingInterval());
        std::cout << "tint set, time goes from: " << getCGTimeMin() << " to: " << getTimeMaxStat() << std::endl;
    }

    statFile.open();

    statFile.getFstream() << Points.begin()->write_variable_names() << std::endl;
    statFile.getFstream() << print_CG() << std::endl;
    
    //std::couts variables
    if (verbosity > 1)
        std::cout << std::endl;
    if (verbosity > 0)
        std::cout << std::endl << printStat() << std::endl;
    
}

template<StatType T>
bool StatisticsVector<T>::readNextDataFile(unsigned int format)
{
    static unsigned int count = 0;
    
    if (count)
        for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
        {
            (*it)->setPreviousPosition((*it)->getPosition());
        }
    
    bool ret_val = DPMBase::readNextDataFile(format);
    
    if (!count)
        for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
        {
            (*it)->setPreviousPosition((*it)->getPosition());
        }
    
    count++;
    return ret_val;
}

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

    // write to .stat file
    statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getTime() << std::endl;
    for (unsigned int i = 0; i < Points.size(); i++)
        statFile.getFstream() << Points[i];
    //std::cout << "S" << Points[2] << std::endl;
    if (getDoGradient())
    {
        statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getTime() << std::endl;
        for (unsigned int i = 0; i < dx.size(); i++)
            statFile.getFstream() << dx[i];
        statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getTime() << std::endl;
        for (unsigned int i = 0; i < dy.size(); i++)
            statFile.getFstream() << dy[i];
        statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getTime() << std::endl;
        for (unsigned int i = 0; i < dz.size(); i++)
            statFile.getFstream() << dz[i];
    }
    
}

template<StatType T>
void StatisticsVector<T>::finishStatistics()
{
    if (getDoTimeAverage())
        write_time_average_statistics();
    statFile.close();
}

template<StatType T>
void StatisticsVector<T>::processStatistics(bool usethese)
{
    if (check_current_time_for_statistics() && usethese)
    {
        if (getMirrorAtDomainBoundary() != 0.0)
        {
            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 (getDoGradient())
                    {
                        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 /= mathsFunc::square(Points[i].Density);
        if (getDoTimeAverage())
        {
            for (unsigned int i = 0; i < timeAverage.size(); i++)
            {
                timeAverage[i] += Points[i];
                if (getDoVariance())
                    timeVariance[i] += Points[i].getSquared();
                if (getDoGradient())
                {
                    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 (getDoVariance())
            {
                timeVariance[i].firstTimeAverage(nTimeAverage);
                timeVariance[i] -= timeAverage[i].getSquared();
            }
            if (getDoGradient())
            {
                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 (getDoVariance())
            {
                timeVariance[i] /= nTimeAverage;
                timeVariance[i] -= timeAverage[i].getSquared();
            }
            if (getDoGradient())
            {
                dxTimeAverage[i] /= nTimeAverage;
                dyTimeAverage[i] /= nTimeAverage;
                dzTimeAverage[i] /= nTimeAverage;
            }
        }
    }
    
    // write average to .stat file
    statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getCGTimeMin() << " " << getTime() << std::endl;
    for (unsigned int i = 0; i < timeAverage.size(); i++)
        statFile.getFstream() << timeAverage[i];
    // write to std::cout
    StatisticsPoint<T> avg = average(timeAverage);
    std::cout << "Averages: " << avg.print() << std::endl;
    
    if (getDoVariance())
    {
        // write variance to .stat file
        statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getCGTimeMin() << " " << getTime() << " " << std::endl;
        for (unsigned int i = 0; i < timeVariance.size(); i++)
            statFile.getFstream() << timeVariance[i];
        //StatisticsPoint<T> var = average(timeVariance);
    } //end if (getDoVariance)
    
    if (getDoGradient())
    {
        // write variance to .stat file
        statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getCGTimeMin() << " " << getTime() << " " << std::endl;
        for (unsigned int i = 0; i < timeAverage.size(); i++)
            statFile.getFstream() << dxTimeAverage[i];
        statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getCGTimeMin() << " " << getTime() << " " << std::endl;
        for (unsigned int i = 0; i < timeAverage.size(); i++)
            statFile.getFstream() << dyTimeAverage[i];
        statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getCGTimeMin() << " " << getTime() << " " << std::endl;
        for (unsigned int i = 0; i < timeAverage.size(); i++)
            statFile.getFstream() << dzTimeAverage[i];
    }
    
    setCGTimeMin(getTime());
    
}
template<StatType T>
void StatisticsVector<T>::statistics_from_fstat_and_data()
{
    std::cout << "set " << " positions" << std::endl;
    dataFile.setCounter(0);
    fStatFile.setCounter(0);

    //This opens the file the data will be recalled from
    if (dataFile.getFileType() == FileType::ONE_FILE)
        dataFile.open(std::fstream::in);
    if (fStatFile.getFileType() == FileType::ONE_FILE)
        fStatFile.open(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 (readNextDataFile(format))
            break;
        fStatFile.openNextFile(std::fstream::in);
    }

    //set tminstat to smallest time evaluated
    if (getTime() >= getCGTimeMin())
        setCGTimeMin(getTime() - getTimeStep());

    //This creates the file statistics will be saved to
    //getStatFile.open(std::fstream::out);
    
    if (verbosity > 1)
    {
        std::cout << "Input from " << dataFile.getName() << std::endl;
        std::cout << "Input from " << fStatFile.getName() << std::endl;
        std::cout << "Output to " << statFile.getName() << std::endl << std::endl;
    }
    
    initialiseStatistics();
    
    //std::couts collision time
//    if (verbosity > 1)
//    {
//        //getLargestParticle()->computeMass(getSpecies()); //This should already have been done in readNextdataFile;
//        Mdouble mass = particleHandler.getLargestParticle()->getMass();
//        std::cout << "Collision Time " << speciesHandler.getObject(0)->getCollisionTime(mass)
//            << " and restitution coefficient " << speciesHandler.getObject(0)->getRestitutionCoefficient(mass)
//                << " for mass " << mass << std::endl
//                << std::endl;
//    }
    //std::couts particles
    if (verbosity > 2)
    {
        DPMBase::write(std::cout, true);
        std::cout << std::endl;
    }
    else if (verbosity)
    {
        //DPMBase::write(std::cout, false);
        //std::cout << std::endl;
    }


    //this increases the file counter
    if ((fStatFile.getFileType() == FileType::MULTIPLE_FILES || fStatFile.getFileType() == FileType::MULTIPLE_FILES_PADDED) && getTime() < getCGTimeMin())
    {
        if (verbosity > 1)
            findNextExistingDataFile(getCGTimeMin());
        else
            findNextExistingDataFile(getCGTimeMin(), false);
        for (unsigned int i = 1; i < dataFile.getCounter(); i++)
            jump_fstat();

    }

    //Output statistics for each time step
    std::cout << "Start statistics" << std::endl;
    do
    {
        if (getTime() > getTimeMaxStat())
        {
            std::cout
                << "reached end t=" << std::setprecision(9) << getTime()
                << " tmaxStat=" << std::setprecision(9) << getTimeMaxStat()
                << std::endl;
            break;
        }
        else if (getTime() < getCGTimeMin())
        {
            if (verbosity > 1)
                std::cout << "Jumping statistics t=" << getTime() << " tminStat()=" << getCGTimeMin() << std::endl;
            jump_fstat();
        }
        else
        {
            if (verbosity)
                std::cout << "Get statistics for t=" << std::setprecision(9) << getTime() << std::endl;
            if (verbosity > 1)
                std::cout << "#particles=" << particleHandler.getNumberOfObjects() << std::endl;
            
            if (!walls)
                wallHandler.clear();
            if (!periodicWalls)
                boundaryHandler.clear();
            
            outputStatistics();
            gather_force_statistics_from_fstat_and_data();
            processStatistics(true);
            
        }
        //if step size>1, skip a number of steps.
        if (dataFile.getFileType() == FileType::ONE_FILE)
            for (unsigned int i = 1; i < getStepSize(); i++)
            {
                //you only get here if you change the step size and you use a single data file
                readNextDataFile(format);
                jump_fstat();
            }
    } while (readNextDataFile(format));
    //set tmax to largest time evaluated
    if (getTime() < getTimeMaxStat())
        setTimeMaxStat(getTime());
    
    finishStatistics();
    dataFile.close();
    fStatFile.close();
}

template<StatType T>
void StatisticsVector<T>::auto_setdim()
{
    if (particleHandler.getNumberOfObjects() < 1)
        return;
    std::vector<BaseParticle*>::iterator it = particleHandler.begin();
    setXMin((*it)->getPosition().X - (*it)->getRadius());
    setYMin((*it)->getPosition().Y - (*it)->getRadius());
    setZMin((*it)->getPosition().Z - (*it)->getRadius());
    setXMax((*it)->getPosition().X + (*it)->getRadius());
    setYMax((*it)->getPosition().Y + (*it)->getRadius());
    setZMax((*it)->getPosition().Z + (*it)->getRadius());
    for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
    {
        setXMin(std::min(getXMin(), (*it)->getPosition().X - (*it)->getRadius()));
        setYMin(std::min(getYMin(), (*it)->getPosition().Y - (*it)->getRadius()));
        setZMin(std::min(getZMin(), (*it)->getPosition().Z - (*it)->getRadius()));
        setXMax(std::max(getXMax(), (*it)->getPosition().X + (*it)->getRadius()));
        setYMax(std::max(getYMax(), (*it)->getPosition().Y + (*it)->getRadius()));
        setZMax(std::max(getZMax(), (*it)->getPosition().Z + (*it)->getRadius()));
    } //end for all particles
}

template<StatType T>
bool StatisticsVector<T>::read_next_from_p3p_file()
{
    //uncomment if you want data file output
    //MD::outputXBallsData();
    
    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);
    if (verbosity>2)
        std::cout << dummyStr << std::endl;
    //Reading second line
    p3p_file >> dummy;
    setTime(dummy);
    p3p_file >> N;
    std::cout << "t= " << getTime() << ": loading " << N << " particles ..." << std::endl;
    getline(p3p_file, dummyStr);
    //std::cout << dummyStr << std::endl;
    //Ignoring third line
    getline(p3p_file, dummyStr);
    if (verbosity>2)
        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 (particleHandler.getNumberOfObjects() < N)
        while (particleHandler.getNumberOfObjects() < N)
            particleHandler.copyAndAddObject(P0);
    else
        while (particleHandler.getNumberOfObjects() > N)
            particleHandler.removeLastObject();
    
    //Read forth to (N+3)rd line
    Vec3D dummyVec;
    Vec3D position, velocity, angle, angularVelocity;
    int dummyInt;
    for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
    {
        p3p_file >> dummyInt;
        (*it)->setIndex(dummyInt); //ID
        //std::cout << "  ID " << (*it)->getIndex() << std::endl;
        p3p_file >> dummyInt;
        (*it)->setIndSpecies(dummyInt - 1); //Group
        p3p_file >> dummy;
        (*it)->setRadius(pow(dummy / (4. / 3. * constants::pi), 1. / 3.)); //Volume
        p3p_file >> dummy;
        (*it)->setMass(dummy); //Mass
        p3p_file
                >> position
                >> velocity
                >> angularVelocity;
        (*it)->setPosition(position);
        (*it)->setVelocity(velocity);
        (*it)->setAngularVelocity(angularVelocity);
        //clean up tangential springs
    } //end for all particles
      //auto_setSystemDimensions();
    
    getline(p3p_file, dummyStr);
    if (verbosity>2)
        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 = getName().c_str();
    p3p_filename += ".p3p";
    bool opened = helpers::openFile(p3p_file, p3p_filename, std::fstream::in);
    std::string p3c_filename = getName().c_str();
    p3c_filename += ".p3c";
    helpers::openFile(p3c_file, p3c_filename, std::fstream::in);
    std::string p3w_filename = getName().c_str();
    p3w_filename += ".p3w";
    helpers::openFile(p3w_file, p3w_filename, 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 = getName().c_str();
        p3p_filename += ".p4p";
        opened = helpers::openFile(p3p_file, p3p_filename, std::fstream::in);
        p3c_filename = getName().c_str();
        p3c_filename += ".p4c";
        opened = helpers::openFile(p3c_file, p3c_filename, std::fstream::in);
        p3w_filename = getName().c_str();
        p3w_filename += ".p4w";
        opened = helpers::openFile(p3w_file, p3w_filename, 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 (getTime() >= getCGTimeMin())
        setCGTimeMin(getTime() - getTimeStep());
    DPMBase::writeRestartFile();
    
    dataFile.open(std::fstream::out);
    DPMBase::setSystemDimensions(3);
    setParticleDimensions(3);
    DPMBase::outputXBallsData(dataFile.getFstream());
    DPMBase::writeXBallsScript();
    //exit(-1);
    
    //This creates the file statistics will be saved to
    //getStatFile.open(std::fstream::out);
    
    initialiseStatistics();

    //std::couts collision time
    if (verbosity > 1)
    {
        //getLargestParticle()->computeMass(getSpecies()); //This should already have been done in read_next_from_p3p_file();
        Mdouble mass = particleHandler.getLargestParticle()->getMass();
        auto species = dynamic_cast<LinearViscoelasticSpecies*>(speciesHandler.getObject(0));
        if (species!= nullptr)
        {
            std::cout << "Collision Time " << species->getCollisionTime(mass)
                << " and restitution coefficient " << species->getRestitutionCoefficient(mass)
                << " for mass " << mass << std::endl
                << std::endl;
        }
        //std::couts particles
    }
    if (verbosity > 2)
    {
        DPMBase::write(std::cout, true);
        std::cout << std::endl;
    }
    else if (verbosity)
    {
        DPMBase::write(std::cout,false);
        std::cout << std::endl;
    }

    //Output statistics for each time step
    std::cout << "Start statistics" << std::endl;
    do
    {
        if (getTime() > getTimeMaxStat())
        {
            std::cout << "reached end t=" << std::setprecision(9) << getTime() << " tmaxStat=" << std::setprecision(9) << getTimeMaxStat() << std::endl;
            break;
        }
        else if (getTime() >= getCGTimeMin())
        {
            if (verbosity)
                std::cout << "Get statistics for t=" << std::setprecision(9) << getTime() << std::endl;
            if (verbosity > 1)
                std::cout << "#particles=" << particleHandler.getNumberOfObjects() << std::endl;
            
            if (!walls)
                wallHandler.clear();
            if (!periodicWalls)
                boundaryHandler.clear();
            
            outputStatistics();
            gather_force_statistics_from_p3c(version);
            processStatistics(true);
            
        }
        else
        {
            if (verbosity > 1)
                std::cout
                << "Jumping statistics t=" << getTime()
                        << " tminStat()=" << getCGTimeMin() << std::endl;
            jump_p3c();
        }
        //if step size>1, skip a number of steps.
        for (unsigned int i = 1; i < getStepSize(); 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 (getTime() < getTimeMaxStat())
        setTimeMaxStat(getTime());
    
    finishStatistics();
    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;
    setTime(dummy);
    p3c_file >> N;
    std::cout << "ignoring t= " << getTime() << ": 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);
    if (verbosity>2)
        std::cout << dummyStr << std::endl;
    //Reading second line
    p3w_file >> dummy;
    setTime(dummy);
    p3w_file >> N;
    std::cout << "ignoring t= " << getTime() << ": 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);
    if (verbosity>2)
        std::cout << dummyStr << std::endl;
    //Reading second line
    p3c_file >> dummy;
    setTime(dummy);
    p3c_file >> N;
    std::cout << "t= " << getTime() << ": loading " << N << " contacts ..." << std::endl;
    getline(p3c_file, dummyStr);
    //Ignoring third line
    getline(p3c_file, dummyStr);
    if (verbosity>2)
        std::cout << dummyStr << std::endl;
    //Checking for end of file
    //if (p3c_file.eof()||p3c_file.peek()==-1) return false;
    
    //set up index list
    unsigned int Nmax = 0;
    for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
        Nmax = std::max(Nmax, (*it)->getIndex());
    static std::vector<int> index;
    index.resize(Nmax+1);
    {
        int i = 0;
        for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
        {
            index[(*it)->getIndex()] = 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 //p4p
        {
            p3c_file
                    >> index1_p3
                    >> index2_p3
                    >> Contact
                    >> Force;
        }
        index1 = index[index1_p3];
        index2 = index[index2_p3];
        P1 = particleHandler.getObject(index1);
        P2 = particleHandler.getObject(index2);
        P1_P2_normal = P1->getPosition() - P2->getPosition();
        // deal with periodic boundaries; correct if the distance between two particles is bigger than half the domain size
        // (only works if periodic walls are >2d apart)
        std::cout << "#forces="  << std::setprecision(9) << P1_P2_normal.Y << "|" << getYMax()-getYMin() << std::endl;
        if (2.0*P1_P2_normal.X > getXMax()-getXMin())
            P1_P2_normal.X -=     getXMax()-getXMin();
        else if (2.0*P1_P2_normal.X < getXMin()-getXMax())
            P1_P2_normal.X +=     getXMax()-getXMin();

        else if (2.0*P1_P2_normal.Y > getYMax()-getYMin())
            P1_P2_normal.Y -=     getYMax()-getYMin();
        else if (2.0*P1_P2_normal.Y < getYMin()-getYMax())
            P1_P2_normal.Y +=     getYMax()-getYMin();

        else if (2.0*P1_P2_normal.Z > getZMax()-getZMin())
            P1_P2_normal.Z -=     getZMax()-getZMin();
        else if (2.0*P1_P2_normal.Z < getZMin()-getZMax())
            P1_P2_normal.Z +=     getZMax()-getZMin();
        std::cout << "#forces="  << std::setprecision(9) << P1_P2_normal.Y << "|" << getYMax()-getYMin() << std::endl;

        dist = Vec3D::getLength(P1_P2_normal);
        P1_P2_normal /= dist;
        delta = P1->getRadius() + P2->getRadius() - dist;
        if (version == 3)
        {
            Contact = P2->getPosition() + P1_P2_normal * (P2->getRadius() - delta / 2);
        }
        fdotn = -Vec3D::dot(Force, P1_P2_normal) / Vec3D::getLengthSquared(P1_P2_normal);
        P1_P2_tangential = Force + fdotn * P1_P2_normal; //get tangential force
        fdott = Vec3D::getLength(P1_P2_tangential);
        if (fdott == 0)
        { 
            P1_P2_tangential = Vec3D(0, 0, 0);
        }
        else
        {
            P1_P2_tangential /= fdott;
        }
        //Finished reading file
        std::cout << "#forces=" << P1_P2_normal << "|" << fdotn << std::endl;
        gatherContactStatistics(index1, index2, Contact, delta, 0, fdotn, fdott, P1_P2_normal, -P1_P2_tangential);
        gatherContactStatistics(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);
    if (verbosity>2)
        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);
    if (verbosity>2)
        std::cout << dummyStr << std::endl;
    //Reading second line
    p3w_file >> dummy;
    setTime(dummy);
    p3w_file >> N;
    std::cout << "t= " << getTime() << ": loading " << N << " wall contacts ..." << std::endl;
    getline(p3w_file, dummyStr);
    //Ignoring third line
    getline(p3w_file, dummyStr);
    if (verbosity>2)
        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 = particleHandler.getObject(index1);
        if (version == 3)
        {
            //Branch *= -1;
            Contact = P1->getPosition() + Branch;
        }
        else
        {
            Branch = Contact - P1->getPosition();
        }
        dist = Vec3D::getLength(Branch);
        P1_P2_normal = Branch / dist;
        delta = P1->getRadius() - dist; //check
        fdotn = -Vec3D::dot(Force, P1_P2_normal) / Vec3D::getLengthSquared(P1_P2_normal);
        P1_P2_tangential = Force + fdotn * P1_P2_normal; //get tangential force
        fdott = Vec3D::getLength(P1_P2_tangential);
        if (fdott == 0)
        { 
            P1_P2_tangential = Vec3D(0, 0, 0);
        }
        else
        {
            P1_P2_tangential /= fdott;
        }
        //Finished reading fstat file
        gatherContactStatistics(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);
    if (verbosity>2)
        std::cout << dummyStr << std::endl;
}


template<StatType T>
void StatisticsVector<T>::setPositions()
{
    std::cout << "set " << " positions" << std::endl;
    int N = positions_.size();
    if (N!=0)
    {
        std::cout << "setting " << N << " positions" << std::endl;
        Points.resize(N);
        for (int i = 0; i<N; i++) {
            Points[i].setPosition(positions_[i]);
        }
        nx = N; ny=1; nz=1;
    }
    else //set automatically on mesh
    {
        //Dimensions of the domain / n
        Vec3D diff = Vec3D(
            (getXMaxStat() - getXMinStat()) / nx,
            (getYMaxStat() - getYMinStat()) / ny,
            (getZMaxStat() - getZMinStat()) / nz);
        //Center of the domain
        Vec3D avg = 0.5 * Vec3D(
            getXMaxStat() + getXMinStat(),
            getYMaxStat() + getYMinStat(),
            getZMaxStat() + getZMinStat());
        //Left endpoint of the domain
        Vec3D Min = Vec3D(
            getXMinStat(),
            getYMinStat(),
            getZMinStat());

        int num[3] = {0, 0, 0};
        if (getMirrorAtDomainBoundary() != 0.0)
        {
            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] = static_cast<int>(std::floor(std::max(getMirrorAtDomainBoundary(), getCutoff()) / diff.X));
                nx += 2 * num[0];
                Min.X -= num[0] * diff.X;
            }
            if (statType == Y || statType == XY || statType == YZ || statType == XYZ)
            {
                num[1] = static_cast<int>(std::floor(std::max(getMirrorAtDomainBoundary(), getCutoff()) / diff.Y));
                ny += 2 * num[1];
                Min.Y -= num[1] * diff.Y;
            }
            if (statType == Z || statType == XZ || statType == XZ || statType == XYZ)
            {
                num[2] = static_cast<int>(std::floor(std::max(getMirrorAtDomainBoundary(), getCutoff()) / diff.Z));
                nz += 2 * num[2];
                Min.Z -= num[2] * diff.Z;
            }
        }

        //Set position of statistics points and set variables to 0
        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].setPosition(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 (doDoublePoints)
                            Points[n].setPosition(Points[n].getPosition()
                                - Vec3D((i % 2) ? (0.99 * diff.X) : 0, (j % 2) ? (0.99 * diff.Y) : 0, (k % 2) ? (0.99 * diff.Z) : 0));
                        if (getMirrorAtDomainBoundary() != 0.0) 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].setPosition(Points[k].getPosition().getFromCylindricalCoordinates());
            }
        }
    }
    for (int i = 0; i<N; i++)
    {
        Points[i].setCGInverseVolume();
        Points[i].set_zero();
    }
    if (getDoGradient())
    {
        dx.resize(N);
        dy.resize(N);
        dz.resize(N);
        for (int n = 0; n < N; n++)
        {
            dx[n].setPosition(Points[n].getPosition());
            dx[n].mirrorParticle = Points[n].mirrorParticle;
            dx[n].set_zero();
            dy[n].setPosition(Points[n].getPosition());
            dy[n].mirrorParticle = Points[n].mirrorParticle;
            dy[n].set_zero();
            dz[n].setPosition(Points[n].getPosition());
            dz[n].mirrorParticle = Points[n].mirrorParticle;
            dz[n].set_zero();
        }
    }

    //Set TimeAverage the same way
    if (getDoTimeAverage())
    {
        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].setPosition(Points[n].getPosition());
                        timeAverage[n].set_zero();
                        timeAverage[n].mirrorParticle = Points[n].mirrorParticle;
                        n++;
                    }
        }

        if (getDoVariance())
        {
            timeVariance.resize(N);
            for (int n = 0; n < N; n++)
            {
                timeVariance[n].setPosition(timeAverage[n].getPosition());
                timeVariance[n].mirrorParticle = Points[n].mirrorParticle;
                timeVariance[n].set_zero();
            }
        }

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

    } //end if (getDoTimeAverage())
}

template<StatType T>
bool StatisticsVector<T>::loadPositions(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 positions file " << filename << " failed" << std::endl;
        return false;
    }

    double N=0;
    is >> N;
    positions_.resize(N);
    for (int i=0; i<N; i++)
    {
        is >> positions_[i];
    }
    //Close the file
    is.close();
    std::cout << "Loaded " << positions_.size() << " positions from file " << filename << std::endl;
    return true;
}

template<StatType T>
void StatisticsVector<T>::jump_fstat()
{
    if (fStatFile.getFileType() == FileType::MULTIPLE_FILES || fStatFile.getFileType() == FileType::MULTIPLE_FILES_PADDED)
    {
        fStatFile.openNextFile();
        if (verbosity > 1)
            std::cout << "Using fstat file counter: " << fStatFile.getCounter() << std::endl;
    }
    else
    {
        //read in first three lines (should start with '#')
        std::string dummy;
        getline(fStatFile.getFstream(), dummy);
        getline(fStatFile.getFstream(), dummy);
        getline(fStatFile.getFstream(), dummy);
        //read in all lines belonging to this timestep
        while ((fStatFile.getFstream().peek() != -1) && (fStatFile.getFstream().peek() != '#'))
            getline(fStatFile.getFstream(), dummy);
    }
}

template<StatType T>
bool StatisticsVector<T>::satisfiesInclusionCriteria(BaseParticle* P)
{
    return P->getRadius() > rmin
            && P->getRadius() < rmax
            && P->getPosition().Z < hmax
            && (indSpecies<0 || P->getIndSpecies() == indSpecies);
}

template<StatType T>
void StatisticsVector<T>::gatherContactStatistics(unsigned int index1, int index2, Vec3D Contact, Mdouble delta, Mdouble ctheta UNUSED, 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 (particleHandler.getObject(index1)->isFixed())
        {
            int tmp = index1;
            index1 = index2;
            index2 = tmp;
            P1_P2_normal_ *= -1;
        }
        
        if (satisfiesInclusionCriteria(particleHandler.getObject(index1)))
        {
            P1_P2_normal = P1_P2_normal_;
            
            if (getSystemDimensions() == 2)
                Contact.Z = 0.0;
            
            if (index2 < 0)
            { //wall-particle collision
                P1_P2_distance = particleHandler.getObject(index1)->getRadius() - 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 = (particleHandler.getObject(index1)->getVelocity()) / 2;
                P1_P2_VelocityDifference = (particleHandler.getObject(index1)->getVelocity());
                
                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 (particleHandler.getObject(index2)->isFixed() || !satisfiesInclusionCriteria(particleHandler.getObject(index2)))
            { //particle-fixed particle collision (external force acts at contact point)
                if (particleHandler.getObject(index2)->isFixed() && StressTypeForFixedParticles == 3)
                {
                    //infinite extension of stress
                    P1_P2_distance = particleHandler.getObject(index1)->getRadius() + particleHandler.getObject(index2)->getRadius() - 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 || (!particleHandler.getObject(index2)->isFixed() && StressTypeForFixedParticles == 3))
                {
                    //add force from flowing particle to contact point to stress
                    P1_P2_distance = particleHandler.getObject(index1)->getRadius() + particleHandler.getObject(index2)->getRadius() - 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 = particleHandler.getObject(index1)->getRadius() - 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 = particleHandler.getObject(index1)->getRadius() + particleHandler.getObject(index2)->getRadius() - 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 = particleHandler.getObject(index2)->getVelocity() / 2;
                P1_P2_VelocityDifference = -particleHandler.getObject(index2)->getVelocity();
            }
            else
            { //particle-particle collision
            
                if (particleHandler.getObject(index2)->isFixed())
                    std::cout << "ERROR" << std::endl;
                P1_P2_distance = particleHandler.getObject(index1)->getRadius() + particleHandler.getObject(index2)->getRadius() - 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 * particleHandler.getObject(index1)->getRadius() / (particleHandler.getObject(index2)->getRadius() + particleHandler.getObject(index1)->getRadius());
                P1_P2_VelocityAverage = (particleHandler.getObject(index1)->getVelocity() + particleHandler.getObject(index2)->getVelocity()) / 2;
                P1_P2_VelocityDifference = (particleHandler.getObject(index1)->getVelocity() - particleHandler.getObject(index2)->getVelocity());
            }
            //Particles[max(index1,index2)] is chosen so that it is the non-wall, non-fixed particle
            P1_P2_Fabric = MatrixSymmetric3D::selfDyadic(P1_P2_normal) * particleHandler.getObject(index1)->getVolume();
            //fix
            //~ P1_P2_Fabric = SymmetrizedDyadic(P1_P2_normal, P1_P2_normal) * particleHandler.getObject(max(index1,index2))->getVolume(speciesHandler.getObject());
            P1_P2_NormalStress = Matrix3D::dyadic(P1_P2_normal, P1_P2_normal) * (fdotn * norm_dist);
            P1_P2_TangentialStress = Matrix3D::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 = Matrix3D::dyadic(P1_P2_NormalTraction + P1_P2_TangentialTraction, -0.5 * P1_P2_normal * norm_dist);
            P1_P2_Contact = Contact;
            //P1_P2_Potential = (speciesHandler.getObject(0)->getStiffness() * mathsFunc::square(delta) + dynamic_cast<FrictionalSpecies*>(speciesHandler.getObject(0))->getSlidingStiffness() * mathsFunc::square(ctheta)) / 2;
            P1_P2_Potential = 0.0;
            P1_P2_Force = fdotn * P1_P2_normal + fdott * P1_P2_tangential;
            P1_P2_Dissipation = Vec3D::dot(P1_P2_VelocityDifference, P1_P2_Force) / 2;
            P1_P2_CollisionalHeatFlux = P1_P2_normal * (P1_P2_distance / 2 * Vec3D::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].dotNonAveraged(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 || particleHandler.getObject(index2)->isFixed() || !satisfiesInclusionCriteria(particleHandler.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 || particleHandler.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 (fStatFile.getFileType() == FileType::MULTIPLE_FILES || fStatFile.getFileType() == FileType::MULTIPLE_FILES_PADDED)
        fStatFile.openNextFile(std::fstream::in);
    getline(fStatFile.getFstream(), dummy);
    getline(fStatFile.getFstream(), dummy);
    getline(fStatFile.getFstream(), 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 ((fStatFile.getFstream().peek() != -1) && (fStatFile.getFstream().peek() != '#'))
    { //(!fStatFile.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
         */
        fStatFile.getFstream()
                >> time
                >> index1
                >> index2
                >> Contact
                >> delta
                >> ctheta
                >> fdotn
                >> fdott
                >> P1_P2_normal_
                >> P1_P2_tangential;
        fStatFile.getFstream().ignore(256, '\n');
        //Finished reading fstat file
        gatherContactStatistics(index1, index2, Contact, delta, ctheta, fdotn, fdott, P1_P2_normal_, P1_P2_tangential);
        
    }
    if (verbosity > 1)
        std::cout << "#forces=" << counter << ", t=" << time << std::endl;
}

template<StatType T>
void StatisticsVector<T>::evaluate_force_statistics(int wp)
{
    if (periodicWalls)
    {
        for (unsigned int k = wp; k < boundaryHandler.getNumberOfObjects(); k++)
        {
            PeriodicBoundary *b0 = dynamic_cast<PeriodicBoundary *>(boundaryHandler.getObject(k));
            if (b0)
            {
                b0->getDistance(P1);
                b0->shiftPositions(P1, P2);
                evaluate_force_statistics(k + 1);
                b0->shiftPositions(P1, P2);
            }
            AngledPeriodicBoundary* b1 = dynamic_cast<AngledPeriodicBoundary*>(boundaryHandler.getObject(k));
            if (b1 != nullptr)
            {
                b1->distance(P1);
                b1->shiftPositions(P1, P2);
                evaluate_force_statistics(k + 1);
                b1->shiftPositions(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 != 0.0)
        {
            counter++;
            if (!std::isfinite(psi))
            {
                std::cerr << "error: psi =" << psi << " is infinite for "
                    <<   "P1=" << P1
                    << ", P2=" << P2
                    << ", P1_P2_normal=" << P1_P2_normal
                    << ", P1_P2_distance=" << P1_P2_distance
                    << ", t=" << getTime() << std::endl;
                psi = 0;
            }
            Points[i].ContactCoupleStress += (Matrix3D::cross(P1_P2_Contact * psi, P1_P2_ContactCoupleStress) - Matrix3D::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 (getDoGradient())
            {
                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 getDoGradient
        } //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 < boundaryHandler.getNumberOfObjects(); k++)
        {
            PeriodicBoundary* b0 = dynamic_cast<PeriodicBoundary*>(boundaryHandler.getObject(k));
            if (b0)
            {
                b0->getDistance(P1);
                b0->shiftPositions(P1, P2);
                evaluate_force_statistics(k + 1);
                b0->shiftPositions(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 != 0.0)
        {
            Points[i].NormalTraction += P1_P2_NormalTraction * phi;
            Points[i].TangentialTraction += P1_P2_TangentialTraction * phi;
            if (getDoGradient())
            {
                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 getDoGradient
        } //end if phi
    } // end forall Points
    
}

template<StatType T>
void StatisticsVector<T>::outputStatistics()
{
    if (check_current_time_for_statistics())
    {
        for (unsigned int i = 0; i < Points.size(); i++)
            Points[i].setCGInverseVolume();
        
        for (std::vector<BaseParticle*>::iterator P = particleHandler.begin(); P != particleHandler.end(); ++P)
        {
            //unfix particles (turn off to ignore fixed particles)
            if ((!ignoreFixedParticles) && (*P)->isFixed())
            {
                (*P)->unfix();
            }
            //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 < boundaryHandler.getNumberOfObjects(); k++)
        {
            PeriodicBoundary* b = dynamic_cast<PeriodicBoundary*>(boundaryHandler.getObject(k));
            if (b != nullptr)
            {
                Mdouble dist = b->getDistance(P1);
                if (mathsFunc::square(dist - (*P)->getInteractionRadius()) < getCutoff2())
                {
                    b->shiftPosition(*P);
                    evaluate_particle_statistics(P, k + 1);
                    b->shiftPosition(*P);
                }
            }
            AngledPeriodicBoundary* c = dynamic_cast<AngledPeriodicBoundary*>(boundaryHandler.getObject(k));
            if (c != nullptr)
            {
                Mdouble dist = c->distance(P1);
                if (mathsFunc::square(dist - (*P)->getInteractionRadius()) < getCutoff2())
                {
                    c->shiftPosition(*P);
                    evaluate_particle_statistics(P, k + 1);
                    c->shiftPosition(*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)->getPosition());
    
    for (unsigned int i = 0; i < Points.size(); i++)
    {
        phi = Points[i].CG_function((*P)->getPosition());
        if (phi != 0.0)
        {
            if (!std::isfinite(phi))
            {
                std::cout
                << "t =" << DPMBase::getTime()
                        << "error: phi =" << phi
                        << " is infinite for P=" << (*P)->getPosition()
                        << ", Point=" << Points[i].getPosition()
                        << std::endl;
                phi = 0;
            }
            Points[i].Nu += (*P)->getVolume() * phi;
            Points[i].Density += (*P)->getMass() * phi;
            if (VelocityProfile.size())
            {
                Points[i].Momentum += ((*P)->getVelocity() - V) * ((*P)->getMass() * phi);
                Points[i].MomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getVelocity() - V) * ((*P)->getMass() * phi);
            }
            else
            {
                Points[i].Momentum += (*P)->getVelocity() * ((*P)->getMass() * phi);
                Points[i].MomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getVelocity()) * ((*P)->getMass() * phi);
            }
            Points[i].DisplacementMomentum += (*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()) * ((*P)->getMass() * phi);
            Points[i].DisplacementMomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount())) * ((*P)->getMass() * phi);
            Points[i].EnergyFlux += (*P)->getVelocity() * ((*P)->getMass() * (*P)->getVelocity().getLengthSquared() / 2 * phi);
            
            Vec3D LocalAngularMomentum = Vec3D::cross(Points[0].clearAveragedDirections((*P)->getPosition() - Points[i].getPosition()), (*P)->getVelocity()) * (*P)->getMass()
                    + (*P)->getAngularVelocity() * (*P)->getInertia();
            Points[i].LocalAngularMomentum += phi * LocalAngularMomentum;
            Points[i].LocalAngularMomentumFlux += Matrix3D::dyadic(LocalAngularMomentum, (*P)->getVelocity() * phi);
            
            //Note: Displacement is only computed directly if gradients are cmputed
            if (getDoGradient())
            {
                
                Vec3D d_phi = Points[i].CG_gradient((*P)->getPosition(), phi);
                
                if (doDisplacement)
                {
                    for (std::vector<BaseParticle*>::iterator Q = particleHandler.begin(); Q != particleHandler.end(); ++Q)
                    {
                        double phiQ = Points[i].CG_function((*Q)->getPosition()); //Course graining function
                        Points[i].Displacement += MatrixSymmetric3D::symmetrisedDyadic((*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()) - (*Q)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()), d_phi) * ((*P)->getMass() * (*Q)->getMass() * phiQ);
                    }
                    //end:Displacement
                }
                
                dx[i].Nu += (*P)->getVolume() * d_phi.X;
                dx[i].Density += (*P)->getMass() * d_phi.X;
                dx[i].Momentum += (*P)->getVelocity() * ((*P)->getMass() * d_phi.X);
                dx[i].DisplacementMomentum += (*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()) * ((*P)->getMass() * d_phi.X);
                dx[i].MomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getVelocity()) * ((*P)->getMass() * d_phi.X);
                dx[i].DisplacementMomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount())) * ((*P)->getMass() * d_phi.X);
                dx[i].EnergyFlux += (*P)->getVelocity() * ((*P)->getMass() * (*P)->getVelocity().getLengthSquared() / 2 * d_phi.X);
                
                dy[i].Nu += (*P)->getVolume() * d_phi.Y;
                dy[i].Density += (*P)->getMass() * d_phi.Y;
                dy[i].Momentum += (*P)->getVelocity() * ((*P)->getMass() * d_phi.Y);
                dy[i].DisplacementMomentum += (*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()) * ((*P)->getMass() * d_phi.Y);
                dy[i].MomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getVelocity()) * ((*P)->getMass() * d_phi.Y);
                dy[i].DisplacementMomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount())) * ((*P)->getMass() * d_phi.Y);
                dy[i].EnergyFlux += (*P)->getVelocity() * ((*P)->getMass() * (*P)->getVelocity().getLengthSquared() / 2 * d_phi.Y);
                
                dz[i].Nu += (*P)->getVolume() * d_phi.Z;
                dz[i].Density += (*P)->getMass() * d_phi.Z;
                dz[i].Momentum += (*P)->getVelocity() * ((*P)->getMass() * d_phi.Z);
                dz[i].DisplacementMomentum += (*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()) * ((*P)->getMass() * d_phi.Z);
                dz[i].MomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getVelocity()) * ((*P)->getMass() * d_phi.Z);
                dz[i].DisplacementMomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount())) * ((*P)->getMass() * d_phi.Z);
                dz[i].EnergyFlux += (*P)->getVelocity() * ((*P)->getMass() * (*P)->getVelocity().getLengthSquared() / 2 * d_phi.Z);
            } //end if getDoGradient
        } //end if phi
    } // end forall Points
}

//template class StatisticsPoint<XYZ>;

//Templated functions

template<> void StatisticsVector<XY>::setNZ(int new_ UNUSED)
{
    nz = 1;
}
template<> void StatisticsVector<XZ>::setNY(int new_ UNUSED)
{
    ny = 1;
}
template<> void StatisticsVector<YZ>::setNX(int new_ UNUSED)
{
    nx = 1;
}
template<> void StatisticsVector<X>::setNY(int new_ UNUSED)
{
    ny = 1;
}
template<> void StatisticsVector<X>::setNZ(int new_ UNUSED)
{
    nz = 1;
}
template<> void StatisticsVector<Y>::setNX(int new_ UNUSED)
{
    nx = 1;
}
template<> void StatisticsVector<Y>::setNZ(int new_ UNUSED)
{
    nz = 1;
}
template<> void StatisticsVector<Z>::setNX(int new_ UNUSED)
{
    nx = 1;
}
template<> void StatisticsVector<Z>::setNY(int new_ UNUSED)
{
    ny = 1;
}
template<> void StatisticsVector<O>::setNX(int new_ UNUSED)
{
    nx = 1;
}
template<> void StatisticsVector<O>::setNY(int new_ UNUSED)
{
    ny = 1;
}
template<> void StatisticsVector<O>::setNZ(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 - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
            std::max((getXMaxStat() - P2.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P2.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X)),
            std::max((getYMaxStat() - P2.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P2.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y));
}
template<> Mdouble StatisticsVector<RA>::setInfinitelyLongDistance()
{
    return std::min(std::max((getYMaxStat() - P2.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P2.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y),
            std::max((getXMaxStat() - P2.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P2.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X));
}
template<> Mdouble StatisticsVector<RZ>::setInfinitelyLongDistance()
{
    return std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
            std::max((getXMaxStat() - P2.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P2.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X));
}
template<> Mdouble StatisticsVector<AZ>::setInfinitelyLongDistance()
{
    return std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
            std::max((getYMaxStat() - P2.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P2.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y));
}
template<> Mdouble StatisticsVector<R>::setInfinitelyLongDistance()
{
    return std::max((getXMaxStat() - P2.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P2.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X);
}
template<> Mdouble StatisticsVector<A>::setInfinitelyLongDistance()
{
    return std::max((getYMaxStat() - P2.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P2.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y);
}

template<> Mdouble StatisticsVector<XYZ>::setInfinitelyLongDistance()
{
    return std::min(std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
            std::max((getXMaxStat() - P1.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P1.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X)),
            std::max((getYMaxStat() - P1.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P1.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y));
}
template<> Mdouble StatisticsVector<XY>::setInfinitelyLongDistance()
{
    return std::min(std::max((getYMaxStat() - P1.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P1.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y),
            std::max((getXMaxStat() - P1.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P1.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X));
}
template<> Mdouble StatisticsVector<XZ>::setInfinitelyLongDistance()
{
    return std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
            std::max((getXMaxStat() - P1.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P1.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X));
}
template<> Mdouble StatisticsVector<YZ>::setInfinitelyLongDistance()
{
    return std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
            std::max((getYMaxStat() - P1.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P1.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y));
}
template<> Mdouble StatisticsVector<X>::setInfinitelyLongDistance()
{
    return std::max((getXMaxStat() - P1.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P1.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X);
}
template<> Mdouble StatisticsVector<Y>::setInfinitelyLongDistance()
{
    return std::max((getYMaxStat() - P1.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P1.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y);
}
template<> Mdouble StatisticsVector<Z>::setInfinitelyLongDistance()
{
    return fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z);
}
template<> double StatisticsVector<O>::setInfinitelyLongDistance()
{
    return 0;
}

template<StatType T>
void StatisticsVector<T>::writeOutputFiles()
{
    if (statFile.saveCurrentTimestep(getNtimeSteps()))
    {
        outputStatistics();
        DPMBase::gatherContactStatistics();
        processStatistics(true);
    }

    DPMBase::writeOutputFiles();
}