HG-MD  1
Public Member Functions | Private Attributes
StatisticsVector< T > Class Template Reference

This class is used to extract statistical data from MD simulations. More...

#include <StatisticsVector.h>

Inheritance diagram for StatisticsVector< T >:
Inheritance graph
[legend]
Collaboration diagram for StatisticsVector< T >:
Collaboration graph
[legend]

List of all members.

Public Member Functions

void constructor ()
 Sets up all basic things.
void constructor (string name)
 StatisticsVector ()
 Basic constructor only calls constructor()
 StatisticsVector (string name)
 StatisticsVector (StatisticsVector &other)
 Copy constructor.
 StatisticsVector (unsigned int argc, char *argv[])
 Advanced constructor that accepts arguments from the command line.
void readStatArguments (unsigned int argc, char *argv[])
string print ()
 Outputs member variable values to a string.
void set_statType ()
void print_help ()
void set_nx (int new_)
void set_hx (Mdouble hx)
void set_ny (int new_)
void set_hy (Mdouble hy)
void set_nz (int new_)
void set_hz (Mdouble hz)
void set_n (int n)
void set_h (Mdouble h)
int get_nx ()
int get_ny ()
int get_nz ()
void set_tminStat (Mdouble t)
void set_tmaxStat (Mdouble t)
void set_tintStat (Mdouble t)
Mdouble get_tminStat ()
Mdouble get_tmaxStat ()
Mdouble get_tintStat ()
bool check_current_time_for_statistics ()
void set_CG_type (const char *new_)
void set_CG_type (CG new_)
CG get_CG_type ()
void set_n (int nx_, int ny_, int nz_)
void get_n (int &nx_, int &ny_, int &nz_)
void set_w (Mdouble w)
 Set CG variables w2 and CG_invvolume.
void set_w2 (Mdouble new_)
 Set CG variables w2 and CG_invvolume.
Mdouble get_w ()
Mdouble get_w2 ()
Mdouble get_cutoff ()
Mdouble get_cutoff2 ()
string print_CG ()
 Output coarse graining variables.
StatisticsPoint< T > average (vector< StatisticsPoint< T > > &P)
 Output average of statistical variables.
void reset_statistics ()
 Set all statistical variables to zero.
void statistics_from_fstat_and_data ()
 get StatisticsPoint
void set_doTimeAverage (bool new_)
bool get_doTimeAverage ()
void set_StressTypeForFixedParticles (int new_)
int get_StressTypeForFixedParticles ()
void set_infiniteStressForFixedParticles (bool new_)
void set_mirrorAtDomainBoundary (Mdouble new_)
Mdouble get_mirrorAtDomainBoundary ()
void set_doVariance (bool new_)
bool get_doVariance ()
void set_doGradient (bool new_)
bool get_doGradient ()
void set_superexact (bool new_)
bool get_superexact ()
void set_ignoreFixedParticles (bool new_)
bool get_ignoreFixedParticles ()
void verbose ()
void set_verbosity (int new_)
int get_verbosity ()
void set_walls (bool new_)
bool get_walls ()
void set_periodicWalls (bool new_)
bool get_periodicWalls ()
void set_w_over_rmax (Mdouble new_)
Mdouble get_w_over_rmax ()
void set_Positions ()
 Set position of StatisticsPoint points and set variables to 0.
bool read_next_from_data_file (unsigned int format)
void gather_force_statistics_from_fstat_and_data ()
 get force statistics from particle collisions
void evaluate_force_statistics (int wp=0)
 get force statistics
void evaluate_wall_force_statistics (Vec3D P, int wp=0)
 get force statistics (i.e. first particle is a wall particle)
void jump_fstat ()
void initialize_statistics ()
 Initializes statistics, i.e. setting w2, setting the grid and writing the header lines in the .stat file.
void output_statistics ()
 Calculates statistics for Particles (i.e. not collisions)
void gather_statistics_collision (int index1, int index2, Vec3D Contact, Mdouble delta, Mdouble ctheta, Mdouble fdotn, Mdouble fdott, Vec3D P1_P2_normal_, Vec3D P1_P2_tangential)
 Calculates statistics for one collision (can be any kind of collision)
void process_statistics (bool usethese)
 Processes all gathered statistics and resets them afterwards. (Processing means either calculating time averages or writing out statistics)
void finish_statistics ()
 Finish all statistics (i.e. write out final data)
void write_statistics ()
 Writes regular statistics.
void write_time_average_statistics ()
 Writes out time averaged statistics.
void evaluate_particle_statistics (vector< Particle * >::iterator P, int wp=0)
 Calculates statistics for a single Particle.
vector< StatisticsPoint< T > > get_Points ()
Mdouble get_xminStat ()
 Functions to acces and change the domain of statistics.
Mdouble get_yminStat ()
Mdouble get_zminStat ()
Mdouble get_xmaxStat ()
Mdouble get_ymaxStat ()
Mdouble get_zmaxStat ()
void set_xminStat (Mdouble xminStat_)
void set_yminStat (Mdouble yminStat_)
void set_zminStat (Mdouble zminStat_)
void set_xmaxStat (Mdouble xmaxStat_)
void set_ymaxStat (Mdouble ymaxStat_)
void set_zmaxStat (Mdouble zmaxStat_)
int get_nTimeAverage ()
Mdouble setInfinitelyLongDistance ()
void set_Polynomial (vector< Mdouble > new_coefficients, unsigned int new_dim)
void set_Polynomial (Mdouble *new_coefficients, unsigned int num_coeff, unsigned int new_dim)
void set_PolynomialName (const char *new_name)
string get_PolynomialName ()
void set_rmin (Mdouble new_)
void set_rmax (Mdouble new_)
void set_hmax (Mdouble new_)
Mdouble evaluatePolynomial (Mdouble r)
Mdouble evaluateIntegral (Mdouble n1, Mdouble n2, Mdouble t)
template<>
void set_nz (int new_ UNUSED)
template<>
void set_ny (int new_ UNUSED)
template<>
void set_nx (int new_ UNUSED)
template<>
void set_ny (int new_ UNUSED)
template<>
void set_nz (int new_ UNUSED)
template<>
void set_nx (int new_ UNUSED)
template<>
void set_nz (int new_ UNUSED)
template<>
void set_nx (int new_ UNUSED)
template<>
void set_ny (int new_ UNUSED)
template<>
void set_nx (int new_ UNUSED)
template<>
void set_ny (int new_ UNUSED)
template<>
void set_nz (int new_ UNUSED)
template<>
void set_statType ()
template<>
void set_statType ()
template<>
void set_statType ()
template<>
void set_statType ()
template<>
void set_statType ()
template<>
void set_statType ()
template<>
void set_statType ()
template<>
void set_statType ()
template<>
Mdouble setInfinitelyLongDistance ()
template<>
Mdouble setInfinitelyLongDistance ()
template<>
Mdouble setInfinitelyLongDistance ()
template<>
Mdouble setInfinitelyLongDistance ()
template<>
Mdouble setInfinitelyLongDistance ()
template<>
Mdouble setInfinitelyLongDistance ()
template<>
Mdouble setInfinitelyLongDistance ()
template<>
double setInfinitelyLongDistance ()

Private Attributes

StatType statType
 Possible values X,Y,Z,XY,XZ,YZ,XYZ are used to determine if the statistics are averaged; f.e. StatType X is averaged over y and z.
int nx
 Grid size nx,ny,nz (by default the points of evaluation are placed in an grid on the domain [xminStat,xmaxStat]x[yminStat,ymaxStat]x[zminStat,zmaxStat].
int ny
 see nx
int nz
 see nx
Mdouble xminStat
 By default the points of evaluation are placed in an grid on the domain [xminStat,xmaxStat]x[yminStat,ymaxStat]x[zminStat,zmaxStat].
Mdouble xmaxStat
Mdouble yminStat
Mdouble ymaxStat
Mdouble zminStat
Mdouble zmaxStat
int nxMirrored
 extension of grid size from mirrored points
int nyMirrored
int nzMirrored
vector< StatisticsPoint< T > > Points
 A vector that stores the values of the statistical variables at a given position.
vector< StatisticsPoint< T > > dx
 A vector that stores the gradient in x of all statistical variables at a given position.
vector< StatisticsPoint< T > > dy
 A vector that stores the gradient in y of all statistical variables at a given position.
vector< StatisticsPoint< T > > dz
 A vector that stores the gradient in z of all statistical variables at a given position.
vector< StatisticsPoint< T > > timeAverage
 A vector used to sum up all statistical values in Points for time-averaging.
vector< StatisticsPoint< T > > timeVariance
 a vector used to sum up the variance in time of all statistical values
vector< StatisticsPoint< T > > dxTimeAverage
 a vector used to sum up all statistical gradients in dx for time-averaging
vector< StatisticsPoint< T > > dyTimeAverage
 a vector used to sum up all statistical gradients in dy for time-averaging
vector< StatisticsPoint< T > > dzTimeAverage
 a vector used to sum up all statistical gradients in dz for time-averaging
bool doTimeAverage
 Determines if output is averaged over time.
int nTimeAverageReset
 Determines after how many timesteps the time average is reset.
bool doVariance
 Determines if variance is outputted.
bool doGradient
 Determines if gradient is outputted.
int nTimeAverage
 A counter needed to average over time.
CG CG_type
 coarse graining type (Gaussian, Heaviside, Polynomial)
NORMALIZED_POLYNOMIAL< T > Polynomial
 Stores the polynomial, if the cg function is an axisymmetric function polynomial in r.
Mdouble w2
 coarse graining width squared; for HeavisideSphere and Gaussian
Mdouble cutoff
 The distance from the origin at which the cg function vanishes; cutoff=w for HeavisideSphere or Polynomial, 3*w for Gaussian.
Mdouble cutoff2
Mdouble w_over_rmax
 if w is not set manually then w will be set by multiplying this value by the largest particle radius at t=0
Mdouble tminStat
 Statistical output will only be created if t>tminStat.
Mdouble tmaxStat
 Statistical output will only be created if t<tmaxStat.
Mdouble tintStat
 Statistical output will only be created if tmaxStat-tintStat< t< tmaxStat.
Mdouble indSpecies
 defines the species for which statistics are extracted (-1 for all species)
Mdouble rmin
 defines the minimum radius of the particles for which statistics are extracted
Mdouble rmax
 defines the maximum radius of the particles for which statistics are extracted
Mdouble hmax
 defines the maximum height of the particles for which statistics are extracted
bool walls
 Turns off walls before evaluation.
bool periodicWalls
 Turns off periodic walls before evaluation (needed for averaging, because we do not yet check if particle is in domain)
bool ignoreFixedParticles
 Determines if fixed particles contribute to particle statistics (density, ...)
int StressTypeForFixedParticles
 0 no Stress from fixed particles 1 Stress from fixed particles distributed between Contact and flowing Particle COM (default) 2 Stress from fixed particles distributed between fixed and flowing Particle COM 3 Stress from fixed particles extends from flowing particle to infinity
int verbosity
 0 no output 1 basic output (timesteps) 2 full output (number of forces and particles, md and stat parameters)
int format
Mdouble mirrorAtDomainBoundary
bool isMDCLR
bool superexact
 If true, cutoff radius for Gaussian is set to 5*w (from 3*w)
Vec3D P1
 Position of first contact point.
Vec3D P2
 Position of second contact point.
Vec3D P1_P2_normal
 Direction of contact.
Mdouble P1_P2_distance
 Length of contact line.
Matrix3D P1_P2_NormalStress
 Contact stress from normal forces along the line of contact.
Matrix3D P1_P2_TangentialStress
 Contact stress from tangential forces along the line of contact.
Vec3D P1_P2_NormalTraction
 Traction from normal forces at contact of flow with fixed particles or walls.
Vec3D P1_P2_TangentialTraction
 Traction from tangential forces at contact of flow with fixed particles or walls.
MatrixSymmetric3D P1_P2_Fabric
 Fabric.
Vec3D P1_P2_CollisionalHeatFlux
 not yet working
Mdouble P1_P2_Dissipation
 not yet working
Mdouble P1_P2_Potential
 not yet working

Detailed Description

template<StatType T>
class StatisticsVector< T >

This class is used to extract statistical data from MD simulations.

When calling statistics_from_fstat_and_data(), statistical data (such as density, pressure, ...) will be extracted at various points in the domain, aligned in a nx*ny*nz grid.

Set functions can be used to define the dimensions of the grid (default: nx=ny=nz=1) and the type and width of the coarse graining function (default: Gaussian, width w=r_max).

Constructor & Destructor Documentation

template<StatType T>
StatisticsVector< T >::StatisticsVector ( ) [inline]

Basic constructor only calls constructor()

References StatisticsVector< T >::constructor().

{constructor();}
template<StatType T>
StatisticsVector< T >::StatisticsVector ( string  name) [inline]

References StatisticsVector< T >::constructor().

{constructor(name);}
template<StatType T>
StatisticsVector< T >::StatisticsVector ( StatisticsVector< T > &  other)

Copy constructor.

References StatisticsVector< T >::CG_type, StatisticsVector< T >::constructor(), StatisticsVector< T >::cutoff, StatisticsVector< T >::cutoff2, StatisticsVector< T >::doGradient, StatisticsVector< T >::doTimeAverage, StatisticsVector< T >::doVariance, StatisticsVector< T >::dx, StatisticsVector< T >::dxTimeAverage, StatisticsVector< T >::dy, StatisticsVector< T >::dyTimeAverage, StatisticsVector< T >::dz, StatisticsVector< T >::dzTimeAverage, StatisticsVector< T >::format, StatisticsVector< T >::hmax, StatisticsVector< T >::ignoreFixedParticles, StatisticsVector< T >::indSpecies, StatisticsVector< T >::isMDCLR, StatisticsVector< T >::mirrorAtDomainBoundary, StatisticsVector< T >::nTimeAverage, StatisticsVector< T >::nx, StatisticsVector< T >::nxMirrored, StatisticsVector< T >::ny, StatisticsVector< T >::nyMirrored, StatisticsVector< T >::nz, StatisticsVector< T >::nzMirrored, StatisticsVector< T >::periodicWalls, StatisticsVector< T >::Points, StatisticsVector< T >::Polynomial, StatisticsVector< T >::rmax, StatisticsVector< T >::rmin, StatisticsVector< T >::StressTypeForFixedParticles, StatisticsVector< T >::superexact, StatisticsVector< T >::timeAverage, StatisticsVector< T >::timeVariance, StatisticsVector< T >::tintStat, StatisticsVector< T >::tmaxStat, StatisticsVector< T >::tminStat, StatisticsVector< T >::verbosity, StatisticsVector< T >::w2, StatisticsVector< T >::w_over_rmax, StatisticsVector< T >::walls, StatisticsVector< T >::xmaxStat, StatisticsVector< T >::xminStat, StatisticsVector< T >::ymaxStat, StatisticsVector< T >::yminStat, StatisticsVector< T >::zmaxStat, and StatisticsVector< T >::zminStat.

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

        timeAverage = other.timeAverage;
        timeVariance = other.timeVariance;
        dxTimeAverage = other.dxTimeAverage;
        dyTimeAverage = other.dyTimeAverage;
        dzTimeAverage = other.dzTimeAverage;
        doTimeAverage = other.doTimeAverage;
        nTimeAverage = 0;
        for (unsigned int i=0; i<timeAverage.size(); i++) timeAverage[i].set_zero();
        for (unsigned int i=0; i<timeVariance.size(); i++) timeVariance[i].set_zero();
        for (unsigned int i=0; i<dxTimeAverage.size(); i++) dxTimeAverage[i].set_zero();
        for (unsigned int i=0; i<dyTimeAverage.size(); i++) dyTimeAverage[i].set_zero();
        for (unsigned int i=0; i<dzTimeAverage.size(); i++) dzTimeAverage[i].set_zero();
        doVariance = other.doVariance;
        doGradient = other.doGradient;
        CG_type = other.CG_type;
        Polynomial = other.Polynomial;
        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;
        cout << endl << endl << "HELLO" << endl << endl;
}
template<StatType T>
StatisticsVector< T >::StatisticsVector ( unsigned int  argc,
char *  argv[] 
)

Advanced constructor that accepts arguments from the command line. 

{
        //check if filename is given
        if (argc<2) {
                cerr<<"fstatistics.exe filename [-options]"<< endl
                << "Type -help for more information" << endl;
                exit(-1);
        } 
        
        //print help if required
        if (!strcmp(argv[1],"-help")) {
                print_help(); 
        }

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

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

Member Function Documentation

template<StatType T>
StatisticsPoint< T > StatisticsVector< T >::average ( vector< StatisticsPoint< T > > &  P)

Output average of statistical variables.

Output names of statistical variables.

References StatisticsPoint< T >::set_zero().

                                                                              {
        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>
bool StatisticsVector< T >::check_current_time_for_statistics ( ) [inline]

References MD::get_dt(), MD::get_t(), StatisticsVector< T >::get_tmaxStat(), and StatisticsVector< T >::get_tminStat().

{return (get_t()>get_tminStat()&&get_t()<=get_tmaxStat()+get_dt());};
template<StatType T>
void StatisticsVector< T >::constructor ( )

Sets up all basic things.

Reimplemented from MD.

References Gaussian, and StatisticsPoint< T >::set_gb().

Referenced by StatisticsVector< T >::StatisticsVector().

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

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

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

        // set default options
        ignoreFixedParticles = true; //this should be true if you want statistics only of the flow
        StressTypeForFixedParticles = 1;
        verbosity = 1;
        walls = true;
        periodicWalls = true;
        format = 0;
        mirrorAtDomainBoundary = false;
        set_superexact(false);
        
        // time averaging
        set_doTimeAverage(true);
        nTimeAverage = 0;
        nTimeAverageReset = -1;
        //calculate variance
        set_doVariance(false);
        //calculate gradient
        set_doGradient(false);

        // additional stuff
        set_options_stat(1);
}
template<StatType T>
void StatisticsVector< T >::constructor ( string  name) 
                                                 { 
        //call default constructor
        constructor();
        
        // set name
        set_name(name.c_str()); 

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

get force statistics

References Vec3D::X, Vec3D::Y, and Vec3D::Z.

{
        if(periodicWalls)
                for (int k=wp; k<get_NWallPeriodic(); k++)
                {
                        get_WallsPeriodic()[k].distance(P1);
                        get_WallsPeriodic()[k].shift_positions(P1,P2);
                        evaluate_force_statistics(k+1);
                        get_WallsPeriodic()[k].shift_positions(P1,P2);
                }
        //evaluate fields that only depend on CParticle parameters
        Mdouble psi; //Course graining integral
        for (unsigned int i=0; i<Points.size(); i++) {
                psi = Points[i].CG_integral(P1, P2, P1_P2_normal, P1_P2_distance);
                if (psi) {
                        Points[i].NormalStress += P1_P2_NormalStress * psi;
                        Points[i].TangentialStress += P1_P2_TangentialStress * psi;
                        Points[i].Fabric += P1_P2_Fabric * psi;
                        Points[i].CollisionalHeatFlux += P1_P2_CollisionalHeatFlux * psi;
                        Points[i].Dissipation += P1_P2_Dissipation * psi;
                        Points[i].Potential += P1_P2_Potential * psi;
                        if (get_doGradient()) {
                                Vec3D d_psi = Points[i].CG_integral_gradient(P1, P2, P1_P2_normal);

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

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

                                dz[i].NormalStress += P1_P2_NormalStress * d_psi.Z;
                                dz[i].TangentialStress += P1_P2_TangentialStress * d_psi.Z;
                                dz[i].Fabric += P1_P2_Fabric * d_psi.Z;
                                dz[i].CollisionalHeatFlux += P1_P2_CollisionalHeatFlux * d_psi.Z;
                                dz[i].Dissipation += P1_P2_Dissipation * d_psi.Z;
                                dz[i].Potential += P1_P2_Potential * d_psi.Z;
                        } //end if get_doGradient
                } //end if phi
        } // end forall Points
}
template<StatType T>
void StatisticsVector< T >::evaluate_particle_statistics ( vector< Particle * >::iterator  P,
int  wp = 0 
)

Calculates statistics for a single Particle.

References sqr, SymmetrizedDyadic(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

                                                                                          {
        if(periodicWalls)
                for (int k=wp; k<get_NWallPeriodic(); k++)
                {
                        Mdouble dist = get_WallsPeriodic()[k].distance((*P)->get_Position());
                        if (sqr(dist-(*P)->get_Radius())<9.0*get_w2())
                        {
                                get_WallsPeriodic()[k].shift_position((*P)->get_Position());
                                evaluate_particle_statistics(P, k+1);
                                get_WallsPeriodic()[k].shift_position((*P)->get_Position());
                        }
                }

        //evaluate fields that only depend on CParticle parameters
        Mdouble phi; //Course graining function

        for (unsigned int i=0; i<Points.size(); i++) {
                phi = Points[i].CG_function((*P)->get_Position());
                if (phi) {
                        Points[i].Nu += (*P)->get_Volume(get_Species()) * phi;
                        Points[i].Density += (*P)->get_Mass() * phi;
                        Points[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * phi);
                        Points[i].DisplacementMomentum += (*P)->get_Displacement() * ((*P)->get_Mass() * phi);
                        Points[i].MomentumFlux += MatrixSymmetric3D(
                                (*P)->get_Velocity().X*(*P)->get_Velocity().X, 
                                (*P)->get_Velocity().X*(*P)->get_Velocity().Y, 
                                (*P)->get_Velocity().X*(*P)->get_Velocity().Z, 
                                (*P)->get_Velocity().Y*(*P)->get_Velocity().Y, 
                                (*P)->get_Velocity().Y*(*P)->get_Velocity().Z, 
                                (*P)->get_Velocity().Z*(*P)->get_Velocity().Z) * ((*P)->get_Mass() * phi);
                        Points[i].DisplacementMomentumFlux += MatrixSymmetric3D(
                                (*P)->get_Displacement().X*(*P)->get_Displacement().X, 
                                (*P)->get_Displacement().X*(*P)->get_Displacement().Y, 
                                (*P)->get_Displacement().X*(*P)->get_Displacement().Z, 
                                (*P)->get_Displacement().Y*(*P)->get_Displacement().Y, 
                                (*P)->get_Displacement().Y*(*P)->get_Displacement().Z, 
                                (*P)->get_Displacement().Z*(*P)->get_Displacement().Z) * ((*P)->get_Mass() * phi);
                        Points[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * phi);

                        //Note: Displacement is only computed directly if gradients are cmputed
                        if (get_doGradient()) {
                                
                                //begin: Linear displacement, \f$1/(2 \rho^2) \sum_{pq} m_p m_q \phi_q *(v_{pqa} \partial_b \phi_p + v_{pqb} \partial_a \phi_p) \f$     
                                Vec3D d_phi = Points[i].CG_gradient((*P)->get_Position(), phi);
                                
                                for (vector<Particle*>::iterator Q = get_ParticleHandler().begin(); Q!=get_ParticleHandler().end(); ++Q) {
                                        double phiQ = Points[i].CG_function((*Q)->get_Position()); //Course graining function
                                        Points[i].Displacement += SymmetrizedDyadic((*P)->get_Displacement() - (*Q)->get_Displacement(), d_phi) 
                                         * ((*P)->get_Mass() * (*Q)->get_Mass() * phiQ);
                                }
                                //end:Displacement

                                dx[i].Nu += (*P)->get_Volume(get_Species()) * d_phi.X;
                                dx[i].Density += (*P)->get_Mass() * d_phi.X;
                                dx[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * d_phi.X);
                                dx[i].DisplacementMomentum += (*P)->get_Displacement() * ((*P)->get_Mass() * d_phi.X);
                                dx[i].MomentumFlux += MatrixSymmetric3D(
                                        (*P)->get_Velocity().X*(*P)->get_Velocity().X, 
                                        (*P)->get_Velocity().X*(*P)->get_Velocity().Y, 
                                        (*P)->get_Velocity().X*(*P)->get_Velocity().Z, 
                                        (*P)->get_Velocity().Y*(*P)->get_Velocity().Y, 
                                        (*P)->get_Velocity().Y*(*P)->get_Velocity().Z, 
                                        (*P)->get_Velocity().Z*(*P)->get_Velocity().Z) * ((*P)->get_Mass() * d_phi.X);
                                dx[i].DisplacementMomentumFlux += MatrixSymmetric3D(
                                        (*P)->get_Displacement().X*(*P)->get_Displacement().X, 
                                        (*P)->get_Displacement().X*(*P)->get_Displacement().Y, 
                                        (*P)->get_Displacement().X*(*P)->get_Displacement().Z, 
                                        (*P)->get_Displacement().Y*(*P)->get_Displacement().Y, 
                                        (*P)->get_Displacement().Y*(*P)->get_Displacement().Z, 
                                        (*P)->get_Displacement().Z*(*P)->get_Displacement().Z) * ((*P)->get_Mass() * d_phi.X);
                                dx[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * d_phi.X);

                                dy[i].Nu += (*P)->get_Volume(get_Species()) * d_phi.Y;
                                dy[i].Density += (*P)->get_Mass() * d_phi.Y;
                                dy[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * d_phi.Y);
                                dy[i].DisplacementMomentum += (*P)->get_Displacement() * ((*P)->get_Mass() * d_phi.Y);
                                dy[i].MomentumFlux += MatrixSymmetric3D(
                                        (*P)->get_Velocity().X*(*P)->get_Velocity().X, 
                                        (*P)->get_Velocity().X*(*P)->get_Velocity().Y, 
                                        (*P)->get_Velocity().X*(*P)->get_Velocity().Z, 
                                        (*P)->get_Velocity().Y*(*P)->get_Velocity().Y, 
                                        (*P)->get_Velocity().Y*(*P)->get_Velocity().Z, 
                                        (*P)->get_Velocity().Z*(*P)->get_Velocity().Z) * ((*P)->get_Mass() * d_phi.Y);
                                dy[i].DisplacementMomentumFlux += MatrixSymmetric3D(
                                        (*P)->get_Displacement().X*(*P)->get_Displacement().X, 
                                        (*P)->get_Displacement().X*(*P)->get_Displacement().Y, 
                                        (*P)->get_Displacement().X*(*P)->get_Displacement().Z, 
                                        (*P)->get_Displacement().Y*(*P)->get_Displacement().Y, 
                                        (*P)->get_Displacement().Y*(*P)->get_Displacement().Z, 
                                        (*P)->get_Displacement().Z*(*P)->get_Displacement().Z) * ((*P)->get_Mass() * d_phi.Y);
                                dy[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * d_phi.Y);

                                dz[i].Nu += (*P)->get_Volume(get_Species()) * d_phi.Z;
                                dz[i].Density += (*P)->get_Mass() * d_phi.Z;
                                dz[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * d_phi.Z);
                                dz[i].DisplacementMomentum += (*P)->get_Displacement() * ((*P)->get_Mass() * d_phi.Z);
                                dz[i].MomentumFlux += MatrixSymmetric3D(
                                        (*P)->get_Velocity().X*(*P)->get_Velocity().X, 
                                        (*P)->get_Velocity().X*(*P)->get_Velocity().Y, 
                                        (*P)->get_Velocity().X*(*P)->get_Velocity().Z, 
                                        (*P)->get_Velocity().Y*(*P)->get_Velocity().Y, 
                                        (*P)->get_Velocity().Y*(*P)->get_Velocity().Z, 
                                        (*P)->get_Velocity().Z*(*P)->get_Velocity().Z) * ((*P)->get_Mass() * d_phi.Z);
                                dz[i].DisplacementMomentumFlux += MatrixSymmetric3D(
                                        (*P)->get_Displacement().X*(*P)->get_Displacement().X, 
                                        (*P)->get_Displacement().X*(*P)->get_Displacement().Y, 
                                        (*P)->get_Displacement().X*(*P)->get_Displacement().Z, 
                                        (*P)->get_Displacement().Y*(*P)->get_Displacement().Y, 
                                        (*P)->get_Displacement().Y*(*P)->get_Displacement().Z, 
                                        (*P)->get_Displacement().Z*(*P)->get_Displacement().Z) * ((*P)->get_Mass() * d_phi.Z);
                                dz[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * d_phi.Z);
                        } //end if get_doGradient
                } //end if phi
        } // end forall Points
}
template<StatType T>
void StatisticsVector< T >::evaluate_wall_force_statistics ( Vec3D  P,
int  wp = 0 
)

get force statistics (i.e. first particle is a wall particle)

References Vec3D::X, Vec3D::Y, and Vec3D::Z.

{
        if(periodicWalls)
                for (int k=wp; k<get_NWallPeriodic(); k++)
                {
                        get_WallsPeriodic()[k].distance(P1);
                        get_WallsPeriodic()[k].shift_positions(P1,P2);
                        evaluate_force_statistics(k+1);
                        get_WallsPeriodic()[k].shift_positions(P1,P2);
                }

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

}
template<StatType T>
Mdouble StatisticsVector< T >::evaluateIntegral ( Mdouble  n1,
Mdouble  n2,
Mdouble  t 
) [inline]

References StatisticsVector< T >::Polynomial.

                                                                    {
                return Polynomial.evaluateIntegral(n1,n2,t);
        }
template<StatType T>
Mdouble StatisticsVector< T >::evaluatePolynomial ( Mdouble  r) [inline]

References StatisticsVector< T >::Polynomial.

                                              {
                return Polynomial.evaluate(r);
        }
template<StatType T>
void StatisticsVector< T >::finish_statistics ( ) [virtual]

Finish all statistics (i.e. write out final data)

Reimplemented from MD.

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

get force statistics from particle collisions 

{
        //read in first three lines (should start with '#')
        string dummy;
        if (get_options_fstat()==2) increase_counter_fstat(fstream::in);
        getline (get_fstat_file(), dummy);
        getline (get_fstat_file(), dummy);
        getline (get_fstat_file(), dummy);
                
                        
        static Mdouble time;
        static int index1, index2;
        static Vec3D Contact;
        static Mdouble delta, ctheta, fdotn, fdott;
        static Vec3D P1_P2_normal_, P1_P2_tangential;
        static Particle PI;
        static Particle PJ;

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

        }
        if (verbosity>1) cout << "#forces="<< counter << endl;
}
template<StatType T>
void StatisticsVector< T >::gather_statistics_collision ( int index1  ,
int index2  ,
Vec3D Contact  ,
Mdouble delta  ,
Mdouble ctheta  ,
Mdouble fdotn  ,
Mdouble fdott  ,
Vec3D P1_P2_normal_  ,
Vec3D P1_P2_tangential   
) [virtual]

Calculates statistics for one collision (can be any kind of collision)

Todo:
: this is to make the stresses infinitely (i.e. 300 units) long; has to be improved

todo{How does this case need to be treated??}

Todo:
this is to make the stresses infinitely (i.e. 300 units) long
Todo:
{this only works for StatType Z}

this is because wall-particle collisions appear only once in fstat

Todo:
{I, Thomas, removed the effect on different particle sizes bc it gave wrong results in the stress balance; why is this here?}
Todo:
{fabric should be divided by N}
Todo:
{Fabric definition only works for monodispersed particles, i.e. C=(F_xx+F_yy+F_zz)/nu only for monodispersed particles}

Note P1_P2 distance and normal are still in 3D, even for averages

Todo:
check this is not killing the Heaviside CG function. If it is donot use the Heaviside function, the simple solution.

todo{Difference here between direct and indirect statistics}

Reimplemented from MD.

References Dyadic(), sqr, and Vec3D::Z.

{
        Vec3D P1_P2_VelocityAverage, P1_P2_VelocityDifference, P1_P2_Force;
        Mdouble norm_dist;
        if(check_current_time_for_statistics())
        {
                if (index2>=0&&get_ParticleHandler().get_Particle(index2)->isFixed()) {
                        //cerr << "Warning" << endl;
                        int tmp = index1;
                        index1 = index2; 
                        index2 = tmp;
                }
                
                if (get_ParticleHandler().get_Particle(index1)->get_Radius()>rmin && get_ParticleHandler().get_Particle(index1)->get_Radius()<rmax)
                {
                        P1_P2_normal=P1_P2_normal_;
                        
                        if (get_dim()==2) Contact.Z = 0.0;
                                
                        if (index2<0)
                        { //wall-particle collision
                                P1_P2_distance = get_ParticleHandler().get_Particle(index1)->get_Radius()-delta;
                                P1 = Contact;
                                P2 = Contact-P1_P2_normal*P1_P2_distance;
                                norm_dist=P1_P2_distance;
                                P1_P2_VelocityAverage = (get_ParticleHandler().get_Particle(index1)->get_Velocity())/2;
                                P1_P2_VelocityDifference = (get_ParticleHandler().get_Particle(index1)->get_Velocity());
                        
                                if (StressTypeForFixedParticles==3) {
                                        P1_P2_distance += 300;
                                        P1 += 300*P1_P2_normal;
                                        //this is because wall-particle collisions appear only once in fstat
                                        norm_dist=P1_P2_distance;
                                }                       
                        }
                        else if (get_ParticleHandler().get_Particle(index1)->isFixed())
                        { //fixed particle-particle collision (external force acts at contact point)
                                //~ P1_P2_distance = Particles[index2].Radius-delta/2;
                                //~ P1 = Contact;
                                //~ P2 = Contact-P1_P2_normal*P1_P2_distance;
                                //~ P1_P2_distance = Particles[index1].Radius+Particles[index2].Radius-delta;
                                //~ P1 = Contact;
                                //~ P2 = Contact-P1_P2_normal*(0.5*P1_P2_distance);
                                //this is to make the the collision count fully                         
                                if (StressTypeForFixedParticles==3) { 
                                        //infinite extension of stress
                                        P1_P2_distance = get_ParticleHandler().get_Particle(index1)->get_Radius()+get_ParticleHandler().get_Particle(index2)->get_Radius()-delta;
                                        P1 = Contact+P1_P2_normal*(0.5*P1_P2_distance);
                                        P2 = Contact-P1_P2_normal*(0.5*P1_P2_distance);
                                        if (P1_P2_normal.Z>0) {
                                                cerr << "Warning: Force normal upwards on base" << 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();
                                        //P1_P2_distance = 1000;
                                        P1 = P2 + P1_P2_normal*P1_P2_distance;
                                        if (P1_P2_normal.Z>0) cout << P1_P2_distance << " z" << P1.Z  << " f" << fdotn << endl;
                                } else if (StressTypeForFixedParticles==1) {
                                        //add force from flowing particle to contact point to stress
                                        P1_P2_distance = get_ParticleHandler().get_Particle(index1)->get_Radius()+get_ParticleHandler().get_Particle(index2)->get_Radius()-delta;
                                        P1 = Contact;
                                        P2 = Contact-P1_P2_normal*(0.5*P1_P2_distance);
                                        P1_P2_distance /= 2.0;
                                        //\todo{Is this right?}
                                        //~ fdotn *=2;
                                        //~ fdott *=2;
                                } else {
                                        //add force from flowing particle to fixed particle to stress
                                        P1_P2_distance = get_ParticleHandler().get_Particle(index1)->get_Radius()+get_ParticleHandler().get_Particle(index2)->get_Radius()-delta;
                                        P1 = Contact+P1_P2_normal*(0.5*P1_P2_distance);
                                        P2 = Contact-P1_P2_normal*(0.5*P1_P2_distance);
                                }
                                norm_dist=P1_P2_distance;//*2.0*Particles[index1].Radius/(Particles[index2].Radius+Particles[index1].Radius);
                                //This is the length of the contact in particle 1.
                                P1_P2_VelocityAverage = get_ParticleHandler().get_Particle(index2)->get_Velocity()/2;
                                P1_P2_VelocityDifference = -get_ParticleHandler().get_Particle(index2)->get_Velocity();
                        }
                        else
                        {//particle-particle collision

                                if (get_ParticleHandler().get_Particle(index2)->isFixed()) cout << "ERROR" << endl;
                                P1_P2_distance = get_ParticleHandler().get_Particle(index1)->get_Radius()+get_ParticleHandler().get_Particle(index2)->get_Radius()-delta;

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

References StatisticsVector< T >::CG_type.

{return CG_type;};
template<StatType T>
Mdouble StatisticsVector< T >::get_cutoff ( ) [inline]

References StatisticsVector< T >::cutoff.

{return cutoff;};
template<StatType T>
Mdouble StatisticsVector< T >::get_cutoff2 ( ) [inline]

References StatisticsVector< T >::cutoff, and sqr.

{return sqr(cutoff);};
template<StatType T>
bool StatisticsVector< T >::get_doGradient ( ) [inline]

References StatisticsVector< T >::doGradient.

{return doGradient;};
template<StatType T>
bool StatisticsVector< T >::get_doTimeAverage ( ) [inline]

References StatisticsVector< T >::doTimeAverage.

{return doTimeAverage;};
template<StatType T>
bool StatisticsVector< T >::get_doVariance ( ) [inline]

References StatisticsVector< T >::doVariance.

{return doVariance;};
template<StatType T>
bool StatisticsVector< T >::get_ignoreFixedParticles ( ) [inline]

References StatisticsVector< T >::ignoreFixedParticles.

{return ignoreFixedParticles;};
template<StatType T>
Mdouble StatisticsVector< T >::get_mirrorAtDomainBoundary ( ) [inline]

References StatisticsVector< T >::mirrorAtDomainBoundary.

{return mirrorAtDomainBoundary;};
template<StatType T>
void StatisticsVector< T >::get_n ( int &  nx_,
int &  ny_,
int &  nz_ 
) [inline]

References StatisticsVector< T >::nx, StatisticsVector< T >::ny, and StatisticsVector< T >::nz.

{nx_=nx; ny_=ny; nz_=nz;};
template<StatType T>
int StatisticsVector< T >::get_nTimeAverage ( ) [inline]

References StatisticsVector< T >::nTimeAverage.

{return nTimeAverage;};
template<StatType T>
int StatisticsVector< T >::get_nx ( ) [inline]

References StatisticsVector< T >::nx.

{return nx;};
template<StatType T>
int StatisticsVector< T >::get_ny ( ) [inline]

References StatisticsVector< T >::ny.

{return ny;};
template<StatType T>
int StatisticsVector< T >::get_nz ( ) [inline]

References StatisticsVector< T >::nz.

{return nz;};
template<StatType T>
bool StatisticsVector< T >::get_periodicWalls ( ) [inline]

References StatisticsVector< T >::periodicWalls.

{return periodicWalls;};
template<StatType T>
vector<StatisticsPoint<T> > StatisticsVector< T >::get_Points ( ) [inline]

References StatisticsVector< T >::Points.

{return Points;};
template<StatType T>
string StatisticsVector< T >::get_PolynomialName ( ) [inline]

References StatisticsVector< T >::Polynomial.

                                    {   
                return Polynomial.get_name();
        }
template<StatType T>
int StatisticsVector< T >::get_StressTypeForFixedParticles ( ) [inline]

References StatisticsVector< T >::StressTypeForFixedParticles.

{return StressTypeForFixedParticles;};
template<StatType T>
bool StatisticsVector< T >::get_superexact ( ) [inline]

References StatisticsVector< T >::superexact.

{return superexact;};
template<StatType T>
Mdouble StatisticsVector< T >::get_tintStat ( ) [inline]

References StatisticsVector< T >::tintStat.

{return tintStat;};
template<StatType T>
Mdouble StatisticsVector< T >::get_tmaxStat ( ) [inline]

References MD::get_tmax(), and StatisticsVector< T >::tmaxStat.

Referenced by StatisticsVector< T >::check_current_time_for_statistics().

{if(isnan(tmaxStat)) return get_tmax(); else return tmaxStat;};
template<StatType T>
Mdouble StatisticsVector< T >::get_tminStat ( ) [inline]

References StatisticsVector< T >::tminStat.

Referenced by StatisticsVector< T >::check_current_time_for_statistics().

{return tminStat;};
template<StatType T>
int StatisticsVector< T >::get_verbosity ( ) [inline]

References StatisticsVector< T >::verbosity.

{return verbosity;};
template<StatType T>
Mdouble StatisticsVector< T >::get_w ( ) [inline]

References StatisticsVector< T >::w2.

{return sqrt(w2);};
template<StatType T>
Mdouble StatisticsVector< T >::get_w2 ( ) [inline]

References StatisticsVector< T >::w2.

{return w2;};
template<StatType T>
Mdouble StatisticsVector< T >::get_w_over_rmax ( ) [inline]

References StatisticsVector< T >::w_over_rmax.

{return w_over_rmax;};
template<StatType T>
bool StatisticsVector< T >::get_walls ( ) [inline]

References StatisticsVector< T >::walls.

{return walls;};
template<StatType T>
Mdouble StatisticsVector< T >::get_xmaxStat ( ) [inline]

References MD::get_xmax(), and StatisticsVector< T >::xmaxStat.

Referenced by StatisticsVector< T >::set_hx().

{if(isnan(xmaxStat)) return get_xmax(); else return xmaxStat;};
template<StatType T>
Mdouble StatisticsVector< T >::get_xminStat ( ) [inline]

Functions to acces and change the domain of statistics.

References MD::get_xmin(), and StatisticsVector< T >::xminStat.

Referenced by StatisticsVector< T >::set_hx().

{if(isnan(xminStat)) return get_xmin(); else return xminStat;};
template<StatType T>
Mdouble StatisticsVector< T >::get_ymaxStat ( ) [inline]

References MD::get_ymax(), and StatisticsVector< T >::ymaxStat.

Referenced by StatisticsVector< T >::set_hy().

{if(isnan(ymaxStat)) return get_ymax(); else return ymaxStat;};
template<StatType T>
Mdouble StatisticsVector< T >::get_yminStat ( ) [inline]

References MD::get_ymin(), and StatisticsVector< T >::yminStat.

Referenced by StatisticsVector< T >::set_hy().

{if(isnan(yminStat)) return get_ymin(); else return yminStat;};
template<StatType T>
Mdouble StatisticsVector< T >::get_zmaxStat ( ) [inline]

References MD::get_zmax(), and StatisticsVector< T >::zmaxStat.

Referenced by StatisticsVector< T >::set_hz().

{if(isnan(zmaxStat)) return get_zmax(); else return zmaxStat;};
template<StatType T>
Mdouble StatisticsVector< T >::get_zminStat ( ) [inline]

References MD::get_zmin(), and StatisticsVector< T >::zminStat.

Referenced by StatisticsVector< T >::set_hz().

{if(isnan(zminStat)) return get_zmin(); else return zminStat;};
template<StatType T>
void StatisticsVector< T >::initialize_statistics ( ) [virtual]

Initializes statistics, i.e. setting w2, setting the grid and writing the header lines in the .stat file.

This creates the file statistics will be saved to

Reimplemented from MD.

References StatisticsVector< T >::set_Positions(), and StatisticsVector< T >::set_w2().

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

        if (!strcmp(get_stat_filename().c_str(),"")) set_stat_filename(); //set ouput filename to default unless otherwise specified
        if(!open_stat_file(fstream::out)) {cerr << "Problem opening stat file aborting" <<endl; exit(-1);};
                
        get_stat_file() << Points.begin()->write_variable_names() << endl;
        get_stat_file() << print_CG() << endl;
}
template<StatType T>
void StatisticsVector< T >::jump_fstat ( ) 
{
        if (get_options_fstat()==2) {
                increase_counter_fstat(fstream::in);
                if (verbosity>1) cout << "Using fstat file counter: "<<get_file_counter() << endl;
        } else {
                //read in first three lines (should start with '#')
                string dummy;
                getline (get_fstat_file(), dummy);
                getline (get_fstat_file(), dummy);
                getline (get_fstat_file(), dummy);
                //read in all lines belonging to this timestep
                while ( (get_fstat_file().peek() != -1) && (get_fstat_file().peek() != '#') )
                        getline (get_fstat_file(), dummy);
        }
}
template<StatType T>
void StatisticsVector< T >::output_statistics ( ) [virtual]

Calculates statistics for Particles (i.e. not collisions)

Because the domain can change in size some stuff has to be updated

Todo:
{Particles shouldn't be unfixed, in case you do live statistics; please correct}
Todo:
We now have an idea of how to deal with fixed particles better, so this can be improved.

Reimplemented from MD.

                                            {
        if(check_current_time_for_statistics())
        {
                for (unsigned int i=0; i<Points.size(); i++)
                        Points[i].set_CG_invvolume();
                        
                for (vector<Particle*>::iterator P = get_ParticleHandler().begin(); P!=get_ParticleHandler().end(); ++P) {
                        //if (!(i++%10000)) cout << "P" << i << "/" << get_N() << endl;
                        //unfix particles (turn off to ignore fixed particles)
                        if ((!ignoreFixedParticles) && (*P)->isFixed()) {
                                (*P)->unfix(get_Species());
                        }
                        //ignore fixed particles

                        if ((*P)->get_Radius()>rmin && (*P)->get_Radius()<rmax && !(*P)->isFixed()) {
                                evaluate_particle_statistics(P);
                        }

                }
        }
}
template<StatType T>
string StatisticsVector< T >::print ( )

Outputs member variable values to a string.

Outputs StatisticsVector.

References Gaussian, HeavisideSphere, and polynomial.

                                  {
        stringstream ss;
        ss      << "Statistical parameters: name=" << get_name() 
                << ",\n nx=" << nx
                << ", ny=" << ny
                << ", nz=" << nz 
                << ", CG_type=";
        if (get_CG_type()==HeavisideSphere) {
                ss << "HeavisideSphere";
        } else if (get_CG_type()==Gaussian) {
                ss << "Gaussian";
        } else if (get_CG_type()==polynomial) {
                ss << get_PolynomialName();
        } else { cerr << "error in CG_function" << endl; exit(-1); }

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

Output coarse graining variables.

Output names of statistical variables. 

                                     {
        stringstream ss;
        ss << "w " << sqrt(get_w2()) 
                << " dim " << get_dim()
                << " domainStat " << get_xminStat() << " " << get_xmaxStat() 
                << " " << get_yminStat() << " " << get_ymaxStat()
                << " " << get_zminStat() << " " << get_zmaxStat()
                << " n " << nx << " " << ny << " " << nz;
        return ss.str();
}
template<StatType T>
void StatisticsVector< T >::print_help ( ) 
{
        cout << "fstatistics.exe filename [-options]: " << endl
                         << "This program evaluates statistical data from .fstat, .data, and. restart files and writes it into a .stat file" << endl << endl;
        cout << "OPTIONS:" << endl << endl;
        cout << "-stattype [X,Y,Z,XY,XZ,YZ,XYZ]: " << endl
                         << "Used to obtain statistics averaged in all coordinate directions but the ones specified; f.e. -stattype Z yields depth-profiles. Default is XYZ." << endl << endl
             << "-CGtype [HeavisideSphere,Gaussian]: " << endl
                         << "Averaging function; default: Gaussian" << endl << endl
             << "-w,-w_over_rmax [Mdouble]: " << endl
                         << "Averaging width, absolute or in multiples of the radius of the largest particle in the restart file; default: w_over_rmax=1" << endl << endl
             << "-n,nx,ny,nz [uint]: " << 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" << endl << endl
             << "-h,hx,hy,hz [Mdouble]: " << endl
                         << "Defines the amount of grid points by setting the mesh size " << endl << endl
             << "-x,y,z [Mdouble Mdouble]: " << 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)" << endl << endl
             << "-indSpecies [int]: " << endl
             << "-rmin,rmax [Mdouble]: " << endl
             << "-hmax [Mdouble]: " << endl
                         << "Allows to only obtain statistics for specific particles" << endl << endl
             << "-periodicwalls [bool]: " << endl
                         << "neglects periodic walls" << endl << endl
             << "-walls [uint]: " << endl
                         << "only takes into account the first n walls" << endl << endl
             << "-verbosity [uint]: " << endl
                         << "amount of screen output (0 minimal, 1 normal, 2 maximal)" << endl << endl
             << "-verbose: " << endl
                         << "identical to -verbosity 2" << endl << endl
             << "-format [uint]: " << endl
                         << "???" << endl << endl
             << "-tmin,tmax,tint [Mdouble]: " << endl
                         << "redefines the time interval [tmin,tmax] for which statistics are evaluated. Default values from restart file; tint sets tmin to tmax-tint." << endl << endl
             << "-stepsize [uint]: " << endl
                         << "to skip some timesteps. Default: 1 " << endl << endl
             << "-statfilename [string]: " << endl
                         << "to change the name of the output file" << endl << endl
             << "-timeaverage [bool]: " << endl
                         << "To average in time (1) or to print statistics for all time steps (0). Default: 1" << endl << endl
             << "-timevariance [bool]: " << endl
                         << "to print the time variance; only for time averaged data. Default: ?" << endl << endl
             << "-gradient: " << endl
                         << "to plot the first derivative of each statistical value" << endl << endl
             << "-ignoreFixedParticles: " << endl
             << "-StressTypeForFixedParticles: " << endl
                 << "-mirrorAtDomainBoundary [Mdouble]" << endl << endl
                         << "" << endl << endl
                         << "OUTPUT:" << endl << endl
                         << "1-3: Coordinates " << endl
                         << "4: Nu " << endl
                         << "5: Density " << endl
                         << "6-8: MomentumX MomentumY MomentumZ " << endl
                         << "9-11: DisplacementMomentumX DisplacementMomentumY DisplacementMomentumZ " << endl
                         << "12-17: MomentumFluxXX MomentumFluxXY MomentumFluxXZ MomentumFluxYY MomentumFluxYZ MomentumFluxZZ" << endl
                         << "18-23: DisplacementMomentumFluxXX DisplacementMomentumFluxXY DisplacementMomentumFluxXZ DisplacementMomentumFluxYY DisplacementMomentumFluxYZ DisplacementMomentumFluxZZ" << endl
                         << "24-26: EnergyFluxX EnergyFluxY EnergyFluxZ " << endl
                         << "27-35: NormalStressXX NormalStressXY NormalStressXZ NormalStressYX NormalStressYY NormalStressYZ NormalStressZX NormalStressZY NormalStressZZ " << endl
                         << "36-44: TangentialStressXX TangentialStressXY TangentialStressXZ TangentialStressYX TangentialStressYY TangentialStressYZ TangentialStressZX TangentialStressZY TangentialStressZZ " << endl
                         << "45-47: NormalTractionX NormalTractionY NormalTractionZ " << endl
                         << "48-50: TangentialTractionX TangentialTractionY TangentialTractionZ " << endl
                         << "51-56: FabricXX FabricXY FabricXZ FabricYY FabricYZ FabricZZ " << endl
                         << "57-59: CollisionalHeatFluxX CollisionalHeatFluxY CollisionalHeatFluxZ " << endl
                         << "60: Dissipation " << endl
                         << "61: Potential " << endl;
                        exit(-1);
}
template<StatType T>
void StatisticsVector< T >::process_statistics ( bool usethese  ) [virtual]

Processes all gathered statistics and resets them afterwards. (Processing means either calculating time averages or writing out statistics)

Reimplemented from MD.

{
        if (check_current_time_for_statistics()&&usethese)
        {
                if (get_mirrorAtDomainBoundary()) {
                        for (unsigned int i=0; i<Points.size(); i++) {
                                if (Points[i].mirrorParticle>=0) {
                                        //add values to the mirrorParticle
                                        Points[Points[i].mirrorParticle]+=Points[i];
                                        Points[i].set_zero();
                                        if (get_doGradient()) {
                                                dx[Points[i].mirrorParticle]+=dx[i];
                                                dy[Points[i].mirrorParticle]+=dy[i];
                                                dz[Points[i].mirrorParticle]+=dz[i];
                                        }
                                }
                        }
                }
                //~ for (unsigned int i=0; i<Points.size(); i++) 
                        //~ Points[i].Displacement /= sqr(Points[i].Density);
                if (get_doTimeAverage()) 
                {
                        for (unsigned int i=0; i<timeAverage.size(); i++)
                        {
                                timeAverage[i] += Points[i];
                                if (get_doVariance()) 
                                        timeVariance[i] += Points[i].getSquared();
                                if (get_doGradient())
                                {
                                        dxTimeAverage[i] += dx[i];
                                        dyTimeAverage[i] += dy[i];
                                        dzTimeAverage[i] += dz[i];
                                }
                        }
                        nTimeAverage++;
                        if (nTimeAverage==nTimeAverageReset) {
                                cout << "write time average" << 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>
bool StatisticsVector< T >::read_next_from_data_file ( unsigned int  format)

Reimplemented from MD.

References MD::read_next_from_data_file().

                                                                       {
        static unsigned int count = 0; 
        Mdouble PrevTime = get_t();
        if (count) for (vector<Particle*>::iterator it = get_ParticleHandler().begin(); it!=get_ParticleHandler().end(); ++it) 
                (*it)->get_Displacement() = (*it)->get_Position();

        bool ret_val = MD::read_next_from_data_file (format);

        if (nTimeAverage) {
                for (vector<Particle*>::iterator it = get_ParticleHandler().begin(); it!=get_ParticleHandler().end(); ++it) { 
                        //correct for periodic boundaries
                        //~ cout << "P" << it->get_Displacement() << "P" << it->get_Position();
                        (*it)->get_Displacement() = ((*it)->get_Position()-(*it)->get_Displacement());
                        if (get_xmax()>get_xmin() && abs((*it)->get_Displacement().X)>.5*(get_xmax()-get_xmin())) {
                                //cout << (*it)->get_Displacement().X << endl;
                                if ((*it)->get_Displacement().X>0) (*it)->get_Displacement().X -= get_xmax()-get_xmin();
                                else (*it)->get_Displacement().X += get_xmax()-get_xmin();
                        }
                        if (get_ymax()>get_ymin() && abs((*it)->get_Displacement().Y)>.5*(get_ymax()-get_ymin())) {
                                //cout << (*it)->get_Displacement().Y << endl;
                                if ((*it)->get_Displacement().Y>0) (*it)->get_Displacement().Y -= get_ymax()-get_ymin();
                                else (*it)->get_Displacement().Y += get_ymax()-get_ymin();
                        }
                        if (get_zmax()>get_zmin() && abs((*it)->get_Displacement().Z)>.5*(get_zmax()-get_zmin())) {
                                //cout << (*it)->get_Displacement().Z << endl;
                                if ((*it)->get_Displacement().Z>0) (*it)->get_Displacement().Z -= get_zmax()-get_zmin();
                                else (*it)->get_Displacement().Z += get_zmax()-get_zmin();
                        }
                        (*it)->get_Displacement() /= (get_t()-PrevTime);
                }
        } else {
                for (vector<Particle*>::iterator it = get_ParticleHandler().begin(); it!=get_ParticleHandler().end(); ++it) 
                        (*it)->get_Displacement().set_zero();
        }
        count++;
        return ret_val;
}
template<StatType T>
void StatisticsVector< T >::readStatArguments ( unsigned int  argc,
char *  argv[] 
)
Todo:
this should also set the ene,restart filename

References Gaussian, and HeavisideSphere.

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

Set all statistical variables to zero. 

                                           {
        for (unsigned int i=0; i<Points.size(); i++) Points[i].set_zero(); 
        if (get_doGradient()) for (unsigned int i=0; i<Points.size(); i++) {
                dx[i].set_zero(); dy[i].set_zero(); dz[i].set_zero(); 
        }
}
template<StatType T>
void StatisticsVector< T >::set_CG_type ( const char *  new_)

References polynomial.

                                                      {
        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 {
                cerr << "Error: unknown CG_type: " << new_ << endl;
                exit(-1);
        }
        set_PolynomialName(new_);
        set_CG_type(polynomial);
}
template<StatType T>
void StatisticsVector< T >::set_CG_type ( CG  new_)

References polynomial.

                                             {
        if (new_==polynomial&&Polynomial.order()<0) {cerr << "Polynomial is not specified" << endl; exit(-1);} 
        CG_type = new_; 
}
template<StatType T>
void StatisticsVector< T >::set_doGradient ( bool  new_) [inline]

References StatisticsVector< T >::doGradient.

{doGradient = new_;};
template<StatType T>
void StatisticsVector< T >::set_doTimeAverage ( bool  new_) [inline]

References StatisticsVector< T >::doTimeAverage.

{doTimeAverage = new_;};
template<StatType T>
void StatisticsVector< T >::set_doVariance ( bool  new_) [inline]

References StatisticsVector< T >::doVariance.

{doVariance = new_;};
template<StatType T>
void StatisticsVector< T >::set_h ( Mdouble  h) [inline]

References StatisticsVector< T >::set_hx(), StatisticsVector< T >::set_hy(), and StatisticsVector< T >::set_hz().

{ set_hx(h); set_hy(h); set_hz(h);};
template<StatType T>
void StatisticsVector< T >::set_hmax ( Mdouble  new_) [inline]

References StatisticsVector< T >::hmax.

{hmax=new_;}
template<StatType T>
void StatisticsVector< T >::set_hx ( Mdouble  hx) [inline]

References StatisticsVector< T >::get_xmaxStat(), StatisticsVector< T >::get_xminStat(), and StatisticsVector< T >::set_nx().

Referenced by StatisticsVector< T >::set_h().

{ set_nx( ceil((get_xmaxStat()-get_xminStat())/hx) ); };
template<StatType T>
void StatisticsVector< T >::set_hy ( Mdouble  hy) [inline]

References StatisticsVector< T >::get_ymaxStat(), StatisticsVector< T >::get_yminStat(), and StatisticsVector< T >::set_ny().

Referenced by StatisticsVector< T >::set_h().

{ set_ny( ceil((get_ymaxStat()-get_yminStat())/hy) ); };
template<StatType T>
void StatisticsVector< T >::set_hz ( Mdouble  hz) [inline]

References StatisticsVector< T >::get_zmaxStat(), StatisticsVector< T >::get_zminStat(), and StatisticsVector< T >::set_nz().

Referenced by StatisticsVector< T >::set_h().

{ set_nz( ceil((get_zmaxStat()-get_zminStat())/hz) ); };
template<StatType T>
void StatisticsVector< T >::set_ignoreFixedParticles ( bool  new_) [inline]

References StatisticsVector< T >::ignoreFixedParticles.

{ignoreFixedParticles = new_;};
template<StatType T>
void StatisticsVector< T >::set_infiniteStressForFixedParticles ( bool  new_) [inline]

References StatisticsVector< T >::StressTypeForFixedParticles.

{StressTypeForFixedParticles = 2+new_;};
template<StatType T>
void StatisticsVector< T >::set_mirrorAtDomainBoundary ( Mdouble  new_) [inline]

References StatisticsVector< T >::mirrorAtDomainBoundary.

{mirrorAtDomainBoundary = new_;};
template<StatType T>
void StatisticsVector< T >::set_n ( int  n) [inline]

References StatisticsVector< T >::set_nx(), StatisticsVector< T >::set_ny(), and StatisticsVector< T >::set_nz().

{ set_nx(n); set_ny(n); set_nz(n); };
template<StatType T>
void StatisticsVector< T >::set_n ( int  nx_,
int  ny_,
int  nz_ 
) [inline]

References StatisticsVector< T >::nx, StatisticsVector< T >::ny, and StatisticsVector< T >::nz.

{nx=nx_; ny=ny_; nz=nz_;};
template<StatType T>
void StatisticsVector< T >::set_nx ( int  new_) [inline]

References StatisticsVector< T >::nx.

Referenced by StatisticsVector< T >::set_hx(), and StatisticsVector< T >::set_n().

{ nx = new_;};
template<>
void StatisticsVector< YZ >::set_nx ( int new_  UNUSED) 
{ nx = 1; }
template<>
void StatisticsVector< Y >::set_nx ( int new_  UNUSED) 
{ nx = 1; }
template<>
void StatisticsVector< Z >::set_nx ( int new_  UNUSED) 
{ nx = 1; }
template<>
void StatisticsVector< O >::set_nx ( int new_  UNUSED) 
{ nx = 1; }
template<StatType T>
void StatisticsVector< T >::set_ny ( int  new_) [inline]

References StatisticsVector< T >::ny.

Referenced by StatisticsVector< T >::set_hy(), and StatisticsVector< T >::set_n().

{ ny = new_;};
template<>
void StatisticsVector< XZ >::set_ny ( int new_  UNUSED) 
{ ny = 1; }
template<>
void StatisticsVector< X >::set_ny ( int new_  UNUSED) 
{ ny = 1; }
template<>
void StatisticsVector< Z >::set_ny ( int new_  UNUSED) 
{ ny = 1; }
template<>
void StatisticsVector< O >::set_ny ( int new_  UNUSED) 
{ ny = 1; }
template<StatType T>
void StatisticsVector< T >::set_nz ( int  new_) [inline]

References MD::get_dim(), and StatisticsVector< T >::nz.

Referenced by StatisticsVector< T >::set_hz(), and StatisticsVector< T >::set_n().

{ nz = new_; if (get_dim()<3) cerr << "Warning in set_nz: dimension less than 3"<<endl; };
template<>
void StatisticsVector< XY >::set_nz ( int new_  UNUSED) 
{ nz = 1; }
template<>
void StatisticsVector< X >::set_nz ( int new_  UNUSED) 
{ nz = 1; }
template<>
void StatisticsVector< Y >::set_nz ( int new_  UNUSED) 
{ nz = 1; }
template<>
void StatisticsVector< O >::set_nz ( int new_  UNUSED) 
{ nz = 1; }
template<StatType T>
void StatisticsVector< T >::set_periodicWalls ( bool  new_) [inline]

References StatisticsVector< T >::periodicWalls.

Referenced by Statistics().

{periodicWalls = new_;};
template<StatType T>
void StatisticsVector< T >::set_Polynomial ( vector< Mdouble >  new_coefficients,
unsigned int  new_dim 
) [inline]

References StatisticsVector< T >::Polynomial.

                                                                                    {   
                Polynomial.set_polynomial(new_coefficients, new_dim);
        }
template<StatType T>
void StatisticsVector< T >::set_Polynomial ( Mdouble *  new_coefficients,
unsigned int  num_coeff,
unsigned int  new_dim 
) [inline]

References StatisticsVector< T >::Polynomial.

                                                                                                     {  
                Polynomial.set_polynomial(new_coefficients, num_coeff, new_dim);
        }
template<StatType T>
void StatisticsVector< T >::set_PolynomialName ( const char *  new_name) [inline]

References StatisticsVector< T >::Polynomial.

                                                      { 
                Polynomial.set_name(new_name);
        }
template<StatType T>
void StatisticsVector< T >::set_Positions ( )

Set position of StatisticsPoint points and set variables to 0.

Set position of statistics points and set variables to 0.

References X, Vec3D::X, XY, XYZ, XZ, Y, Vec3D::Y, YZ, Z, and Vec3D::Z.

Referenced by StatisticsVector< T >::initialize_statistics().

{
        //Dimensions of the domain / n 
        Vec3D diff = Vec3D(
                (get_xmaxStat()-get_xminStat())/nx, 
                (get_ymaxStat()-get_yminStat())/ny,
                (get_zmaxStat()-get_zminStat())/nz);
        //Center of the domain
        Vec3D avg = 0.5 * Vec3D(
                get_xmaxStat()+get_xminStat(), 
                get_ymaxStat()+get_yminStat(), 
                get_zmaxStat()+get_zminStat());
        //Left endpoint of the domain
        Vec3D Min = Vec3D(
                get_xminStat(), 
                get_yminStat(), 
                get_zminStat());

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

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

        //Set TimeAverage the same way
        if (get_doTimeAverage()) {
                timeAverage.resize(N);
                int n = 0;
                for (int i=0; i<max(nx,1); i++)
                        for (int j=0; j<max(ny,1); j++)
                                for (int k=0; k<max(nz,1); k++) {
                                        timeAverage[n].set_Position( Vec3D( nx?(Min.X+diff.X*(i+0.5)):avg.X, 
                                                                                                                                         ny?(Min.Y+diff.Y*(j+0.5)):avg.Y, 
                                                                                                                                         nz?(Min.Z+diff.Z*(k+0.5)):avg.Z ));
                                        timeAverage[n].set_zero();
                                        timeAverage[n].mirrorParticle=Points[n].mirrorParticle;
                                        n++;
                                }

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

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

        } //end if (get_doTimeAverage()) 
}
template<StatType T>
void StatisticsVector< T >::set_rmax ( Mdouble  new_) [inline]

References StatisticsVector< T >::rmax.

{rmax=new_;}
template<StatType T>
void StatisticsVector< T >::set_rmin ( Mdouble  new_) [inline]

References StatisticsVector< T >::rmin.

{rmin=new_;}
template<StatType T>
void StatisticsVector< T >::set_statType ( )
template<>
void StatisticsVector< XYZ >::set_statType ( )

References XYZ.

{ statType=XYZ; }
template<>
void StatisticsVector< XY >::set_statType ( )

References XY.

{ statType=XY; }
template<>
void StatisticsVector< XZ >::set_statType ( )

References XZ.

{ statType=XZ; }
template<>
void StatisticsVector< YZ >::set_statType ( )

References YZ.

{ statType=YZ; }
template<>
void StatisticsVector< X >::set_statType ( )

References X.

{ statType=X; }
template<>
void StatisticsVector< Y >::set_statType ( )

References Y.

{ statType=Y; }
template<>
void StatisticsVector< Z >::set_statType ( )

References Z.

{ statType=Z; }
template<>
void StatisticsVector< O >::set_statType ( )

References O.

{ statType=O; }
template<StatType T>
void StatisticsVector< T >::set_StressTypeForFixedParticles ( int  new_) [inline]

References StatisticsVector< T >::StressTypeForFixedParticles.

{StressTypeForFixedParticles = new_;};
template<StatType T>
void StatisticsVector< T >::set_superexact ( bool  new_) [inline]

References StatisticsVector< T >::superexact.

{superexact = new_;};
template<StatType T>
void StatisticsVector< T >::set_tintStat ( Mdouble  t) [inline]

References MD::t, and StatisticsVector< T >::tintStat.

{tintStat=t;};
template<StatType T>
void StatisticsVector< T >::set_tmaxStat ( Mdouble  t) [inline]

References MD::t, and StatisticsVector< T >::tmaxStat.

{tmaxStat=t;};
template<StatType T>
void StatisticsVector< T >::set_tminStat ( Mdouble  t) [inline]

References MD::t, and StatisticsVector< T >::tminStat.

{tminStat=t;};
template<StatType T>
void StatisticsVector< T >::set_verbosity ( int  new_) [inline]

References StatisticsVector< T >::verbosity.

{verbosity = new_;};
template<StatType T>
void StatisticsVector< T >::set_w ( Mdouble  w) [inline]

Set CG variables w2 and CG_invvolume.

References StatisticsVector< T >::set_w2(), and sqr.

{set_w2(sqr(w));};
template<StatType T>
void StatisticsVector< T >::set_w2 ( Mdouble  new_)

Set CG variables w2 and CG_invvolume.

References Gaussian, HeavisideSphere, polynomial, and sqr.

Referenced by StatisticsVector< T >::initialize_statistics(), and StatisticsVector< T >::set_w().

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

References StatisticsVector< T >::w_over_rmax.

{w_over_rmax = new_;};
template<StatType T>
void StatisticsVector< T >::set_walls ( bool  new_) [inline]

References StatisticsVector< T >::walls.

{walls = new_;};
template<StatType T>
void StatisticsVector< T >::set_xmaxStat ( Mdouble  xmaxStat_) [inline]

References StatisticsVector< T >::xmaxStat.

{xmaxStat=xmaxStat_;};
template<StatType T>
void StatisticsVector< T >::set_xminStat ( Mdouble  xminStat_) [inline]

References StatisticsVector< T >::xminStat.

{xminStat=xminStat_;};
template<StatType T>
void StatisticsVector< T >::set_ymaxStat ( Mdouble  ymaxStat_) [inline]

References StatisticsVector< T >::ymaxStat.

{ymaxStat=ymaxStat_;};
template<StatType T>
void StatisticsVector< T >::set_yminStat ( Mdouble  yminStat_) [inline]

References StatisticsVector< T >::yminStat.

{yminStat=yminStat_;};
template<StatType T>
void StatisticsVector< T >::set_zmaxStat ( Mdouble  zmaxStat_) [inline]

References StatisticsVector< T >::zmaxStat.

{zmaxStat=zmaxStat_;};
template<StatType T>
void StatisticsVector< T >::set_zminStat ( Mdouble  zminStat_) [inline]

References StatisticsVector< T >::zminStat.

{zminStat=zminStat_;};
template<StatType T>
Mdouble StatisticsVector< T >::setInfinitelyLongDistance ( )
template<>
Mdouble StatisticsVector< XYZ >::setInfinitelyLongDistance ( ) 
                                                                    {
        return min(min(max((get_zmax()-P2.Z+5*sqrt(get_w2()))/P1_P2_normal.Z,-(P2.Z-get_zmin()+5*sqrt(get_w2()))/P1_P2_normal.Z),
                                   max((get_xmax()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xmin()+5*sqrt(get_w2()))/P1_P2_normal.X)),
                                   max((get_ymax()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_ymin()+5*sqrt(get_w2()))/P1_P2_normal.Y));
}
template<>
Mdouble StatisticsVector< XY >::setInfinitelyLongDistance ( ) 
                                                                   {
        return min(max((get_ymax()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_ymin()+5*sqrt(get_w2()))/P1_P2_normal.Y),
                           max((get_xmax()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xmin()+5*sqrt(get_w2()))/P1_P2_normal.X));
}
template<>
Mdouble StatisticsVector< XZ >::setInfinitelyLongDistance ( ) 
                                                                   {
        return min(max((get_zmax()-P2.Z+5*sqrt(get_w2()))/P1_P2_normal.Z,-(P2.Z-get_zmin()+5*sqrt(get_w2()))/P1_P2_normal.Z),
                           max((get_xmax()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xmin()+5*sqrt(get_w2()))/P1_P2_normal.X));
}
template<>
Mdouble StatisticsVector< YZ >::setInfinitelyLongDistance ( ) 
                                                                   {
        return min(max((get_zmax()-P2.Z+5*sqrt(get_w2()))/P1_P2_normal.Z,-(P2.Z-get_zmin()+5*sqrt(get_w2()))/P1_P2_normal.Z),
                           max((get_ymax()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_ymin()+5*sqrt(get_w2()))/P1_P2_normal.Y));
}
template<>
Mdouble StatisticsVector< X >::setInfinitelyLongDistance ( ) 
                                                                  {
        return max((get_xmax()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xmin()+5*sqrt(get_w2()))/P1_P2_normal.X);
}
template<>
Mdouble StatisticsVector< Y >::setInfinitelyLongDistance ( ) 
                                                                  {
        return max((get_ymax()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_ymin()+5*sqrt(get_w2()))/P1_P2_normal.Y);
}
template<>
Mdouble StatisticsVector< Z >::setInfinitelyLongDistance ( ) 
                                                                  {
        return max((get_zmax()-P2.Z+5*sqrt(get_w2()))/P1_P2_normal.Z,-(P2.Z-get_zmin()+5*sqrt(get_w2()))/P1_P2_normal.Z);
}
template<>
double StatisticsVector< O >::setInfinitelyLongDistance ( ) 
                                                                 {
        return 0;
}
template<StatType T>
void StatisticsVector< T >::statistics_from_fstat_and_data ( )

get StatisticsPoint

Todo:
{the ignore-periodic-walls should be automated for averaged dimensions}

References MD::print().

Referenced by Statistics().

{
        //This opens the file the data will be recalled from
        set_data_filename();
        open_data_file(fstream::in);

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

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

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

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

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

        finish_statistics();
        get_data_file().close();
        get_fstat_file().close();
}
template<StatType T>
void StatisticsVector< T >::verbose ( ) [inline]

References StatisticsVector< T >::verbosity.

{verbosity = 2;};
template<StatType T>
void StatisticsVector< T >::write_statistics ( )

Writes regular statistics. 

{ 
        cout << "writing statistics for t=" << get_t() << endl;

        // write to .stat file
        get_stat_file() << left << setprecision(3) << setw(6)  << get_t()  << endl;
        for (unsigned int i=0; i<Points.size(); i++)
                get_stat_file() << Points[i];
        if (get_doGradient()) 
        {
                get_stat_file() << left << setprecision(3) << setw(6)  << get_t()  << endl;
                for (unsigned int i=0; i<dx.size(); i++)
                        get_stat_file() << dx[i];
                get_stat_file() << left << setprecision(3) << setw(6)  << get_t()  << endl;
                for (unsigned int i=0; i<dy.size(); i++)
                        get_stat_file() << dy[i];
                get_stat_file() << left << setprecision(3) << setw(6)  << get_t()  << endl;
                for (unsigned int i=0; i<dz.size(); i++)
                        get_stat_file() << dz[i];
        }
        
}
template<StatType T>
void StatisticsVector< T >::write_time_average_statistics ( )

Writes out time averaged statistics.

References StatisticsPoint< T >::print().

{
                if (verbosity) cout << endl << "averaging " << nTimeAverage << " timesteps " << endl;
                if(nTimeAverage==0)
                {
                        fprintf(stderr, "\n\n\tERROR :: Can not do TimeAverage of 0 timesteps\n\n");
                        exit(-1);
                }
                for (unsigned int i=0; i<timeAverage.size(); i++) { 
                        timeAverage[i].timeAverage(nTimeAverage);
                        if (get_doVariance()) {
                                timeVariance[i].timeAverage(nTimeAverage);
                                timeVariance[i] -= timeAverage[i].getSquared();
                        }
                        if (get_doGradient()) {
                                dxTimeAverage[i].timeAverage(nTimeAverage);
                                dyTimeAverage[i].timeAverage(nTimeAverage);
                                dzTimeAverage[i].timeAverage(nTimeAverage);
                        }
                }
                // write average to .stat file
                get_stat_file() << left << setprecision(3) << setw(6) << get_tminStat() << " "  << get_tmaxStat()  << endl;
                for (unsigned int i=0; i<timeAverage.size(); i++)       get_stat_file() << timeAverage[i];
                // write to cout
                StatisticsPoint<T> avg = average(timeAverage);
                cout << "Averages: " << avg.print() << endl;
                
                if (get_doVariance()) {
                        // write variance to .stat file
                        get_stat_file() << left << setprecision(3) << setw(6) << get_tminStat() << " " << get_t() << " " << endl;
                        for (unsigned int i=0; i<timeVariance.size(); i++)
                                get_stat_file() << timeVariance[i];
                        //StatisticsPoint<T> var = average(timeVariance);
                } //end if (get_doVariance)
                
                if (get_doGradient()) {
                        // write variance to .stat file
                        get_stat_file() << left << setprecision(3) << setw(6) <<  get_tminStat() << " " << get_t() << " " << endl;
                        for (unsigned int i=0; i<timeAverage.size(); i++) get_stat_file() << dxTimeAverage[i];
                        get_stat_file() << left << setprecision(3) << setw(6) <<  get_tminStat() << " " << get_t() << " " << endl;
                        for (unsigned int i=0; i<timeAverage.size(); i++) get_stat_file() << dyTimeAverage[i];
                        get_stat_file() << left << setprecision(3) << setw(6) <<  get_tminStat() << " " << get_t() << " " << endl;
                        for (unsigned int i=0; i<timeAverage.size(); i++) get_stat_file() << dzTimeAverage[i];
                }

}

Member Data Documentation

template<StatType T>
CG StatisticsVector< T >::CG_type [private]

coarse graining type (Gaussian, Heaviside, Polynomial)

Referenced by StatisticsVector< T >::get_CG_type(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::cutoff [private]

The distance from the origin at which the cg function vanishes; cutoff=w for HeavisideSphere or Polynomial, 3*w for Gaussian.

Referenced by StatisticsVector< T >::get_cutoff(), StatisticsVector< T >::get_cutoff2(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::cutoff2 [private]

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
bool StatisticsVector< T >::doGradient [private]

Determines if gradient is outputted.

Referenced by StatisticsVector< T >::get_doGradient(), StatisticsVector< T >::set_doGradient(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
bool StatisticsVector< T >::doTimeAverage [private]

Determines if output is averaged over time.

Referenced by StatisticsVector< T >::get_doTimeAverage(), StatisticsVector< T >::set_doTimeAverage(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
bool StatisticsVector< T >::doVariance [private]

Determines if variance is outputted.

Referenced by StatisticsVector< T >::get_doVariance(), StatisticsVector< T >::set_doVariance(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
vector<StatisticsPoint<T> > StatisticsVector< T >::dx [private]

A vector that stores the gradient in x of all statistical variables at a given position.

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
vector<StatisticsPoint<T> > StatisticsVector< T >::dxTimeAverage [private]

a vector used to sum up all statistical gradients in dx for time-averaging

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
vector<StatisticsPoint<T> > StatisticsVector< T >::dy [private]

A vector that stores the gradient in y of all statistical variables at a given position.

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
vector<StatisticsPoint<T> > StatisticsVector< T >::dyTimeAverage [private]

a vector used to sum up all statistical gradients in dy for time-averaging

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
vector<StatisticsPoint<T> > StatisticsVector< T >::dz [private]

A vector that stores the gradient in z of all statistical variables at a given position.

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
vector<StatisticsPoint<T> > StatisticsVector< T >::dzTimeAverage [private]

a vector used to sum up all statistical gradients in dz for time-averaging

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
int StatisticsVector< T >::format [private]

Reimplemented from MD.

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::hmax [private]

defines the maximum height of the particles for which statistics are extracted

Referenced by StatisticsVector< T >::set_hmax(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
bool StatisticsVector< T >::ignoreFixedParticles [private]

Determines if fixed particles contribute to particle statistics (density, ...)

Referenced by StatisticsVector< T >::get_ignoreFixedParticles(), StatisticsVector< T >::set_ignoreFixedParticles(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::indSpecies [private]

defines the species for which statistics are extracted (-1 for all species)

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
bool StatisticsVector< T >::isMDCLR [private]

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::mirrorAtDomainBoundary [private]

Referenced by StatisticsVector< T >::get_mirrorAtDomainBoundary(), StatisticsVector< T >::set_mirrorAtDomainBoundary(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
int StatisticsVector< T >::nTimeAverage [private]

A counter needed to average over time.

Referenced by StatisticsVector< T >::get_nTimeAverage(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
int StatisticsVector< T >::nTimeAverageReset [private]

Determines after how many timesteps the time average is reset.

template<StatType T>
int StatisticsVector< T >::nx [private]

Grid size nx,ny,nz (by default the points of evaluation are placed in an grid on the domain [xminStat,xmaxStat]x[yminStat,ymaxStat]x[zminStat,zmaxStat].

Referenced by StatisticsVector< T >::get_n(), StatisticsVector< T >::get_nx(), StatisticsVector< T >::set_n(), StatisticsVector< T >::set_nx(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
int StatisticsVector< T >::nxMirrored [private]

extension of grid size from mirrored points

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
int StatisticsVector< T >::ny [private]

see nx

Referenced by StatisticsVector< T >::get_n(), StatisticsVector< T >::get_ny(), StatisticsVector< T >::set_n(), StatisticsVector< T >::set_ny(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
int StatisticsVector< T >::nyMirrored [private]

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
int StatisticsVector< T >::nz [private]

see nx

Referenced by StatisticsVector< T >::get_n(), StatisticsVector< T >::get_nz(), StatisticsVector< T >::set_n(), StatisticsVector< T >::set_nz(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
int StatisticsVector< T >::nzMirrored [private]

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
Vec3D StatisticsVector< T >::P1 [private]

Position of first contact point.

template<StatType T>
Vec3D StatisticsVector< T >::P1_P2_CollisionalHeatFlux [private]

not yet working

template<StatType T>
Mdouble StatisticsVector< T >::P1_P2_Dissipation [private]

not yet working

template<StatType T>
Mdouble StatisticsVector< T >::P1_P2_distance [private]

Length of contact line.

template<StatType T>
MatrixSymmetric3D StatisticsVector< T >::P1_P2_Fabric [private]

Fabric.

template<StatType T>
Vec3D StatisticsVector< T >::P1_P2_normal [private]

Direction of contact.

template<StatType T>
Matrix3D StatisticsVector< T >::P1_P2_NormalStress [private]

Contact stress from normal forces along the line of contact.

template<StatType T>
Vec3D StatisticsVector< T >::P1_P2_NormalTraction [private]

Traction from normal forces at contact of flow with fixed particles or walls.

template<StatType T>
Mdouble StatisticsVector< T >::P1_P2_Potential [private]

not yet working

template<StatType T>
Matrix3D StatisticsVector< T >::P1_P2_TangentialStress [private]

Contact stress from tangential forces along the line of contact.

template<StatType T>
Vec3D StatisticsVector< T >::P1_P2_TangentialTraction [private]

Traction from tangential forces at contact of flow with fixed particles or walls.

template<StatType T>
Vec3D StatisticsVector< T >::P2 [private]

Position of second contact point.

template<StatType T>
bool StatisticsVector< T >::periodicWalls [private]

Turns off periodic walls before evaluation (needed for averaging, because we do not yet check if particle is in domain)

Referenced by StatisticsVector< T >::get_periodicWalls(), StatisticsVector< T >::set_periodicWalls(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
vector<StatisticsPoint<T> > StatisticsVector< T >::Points [private]

A vector that stores the values of the statistical variables at a given position.

Referenced by StatisticsVector< T >::get_Points(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
NORMALIZED_POLYNOMIAL<T> StatisticsVector< T >::Polynomial [private]

Stores the polynomial, if the cg function is an axisymmetric function polynomial in r.

Referenced by StatisticsVector< T >::evaluateIntegral(), StatisticsVector< T >::evaluatePolynomial(), StatisticsVector< T >::get_PolynomialName(), StatisticsVector< T >::set_Polynomial(), StatisticsVector< T >::set_PolynomialName(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::rmax [private]

defines the maximum radius of the particles for which statistics are extracted

Referenced by StatisticsVector< T >::set_rmax(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::rmin [private]

defines the minimum radius of the particles for which statistics are extracted

Todo:
Thomas: maybe this fixed condition should be replaced by a condition function, bool include_statistics_if()

Referenced by StatisticsVector< T >::set_rmin(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
StatType StatisticsVector< T >::statType [private]

Possible values X,Y,Z,XY,XZ,YZ,XYZ are used to determine if the statistics are averaged; f.e. StatType X is averaged over y and z.

template<StatType T>
int StatisticsVector< T >::StressTypeForFixedParticles [private]

0 no Stress from fixed particles 1 Stress from fixed particles distributed between Contact and flowing Particle COM (default) 2 Stress from fixed particles distributed between fixed and flowing Particle COM 3 Stress from fixed particles extends from flowing particle to infinity

Referenced by StatisticsVector< T >::get_StressTypeForFixedParticles(), StatisticsVector< T >::set_infiniteStressForFixedParticles(), StatisticsVector< T >::set_StressTypeForFixedParticles(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
bool StatisticsVector< T >::superexact [private]

If true, cutoff radius for Gaussian is set to 5*w (from 3*w)

Referenced by StatisticsVector< T >::get_superexact(), StatisticsVector< T >::set_superexact(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
vector<StatisticsPoint<T> > StatisticsVector< T >::timeAverage [private]

A vector used to sum up all statistical values in Points for time-averaging.

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
vector<StatisticsPoint<T> > StatisticsVector< T >::timeVariance [private]

a vector used to sum up the variance in time of all statistical values

Referenced by StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::tintStat [private]

Statistical output will only be created if tmaxStat-tintStat< t< tmaxStat.

Referenced by StatisticsVector< T >::get_tintStat(), StatisticsVector< T >::set_tintStat(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::tmaxStat [private]

Statistical output will only be created if t<tmaxStat.

Referenced by StatisticsVector< T >::get_tmaxStat(), StatisticsVector< T >::set_tmaxStat(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::tminStat [private]

Statistical output will only be created if t>tminStat.

Referenced by StatisticsVector< T >::get_tminStat(), StatisticsVector< T >::set_tminStat(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
int StatisticsVector< T >::verbosity [private]

0 no output 1 basic output (timesteps) 2 full output (number of forces and particles, md and stat parameters)

Referenced by StatisticsVector< T >::get_verbosity(), StatisticsVector< T >::set_verbosity(), StatisticsVector< T >::StatisticsVector(), and StatisticsVector< T >::verbose().

template<StatType T>
Mdouble StatisticsVector< T >::w2 [private]

coarse graining width squared; for HeavisideSphere and Gaussian

Referenced by StatisticsVector< T >::get_w(), StatisticsVector< T >::get_w2(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::w_over_rmax [private]

if w is not set manually then w will be set by multiplying this value by the largest particle radius at t=0

Referenced by StatisticsVector< T >::get_w_over_rmax(), StatisticsVector< T >::set_w_over_rmax(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
bool StatisticsVector< T >::walls [private]

Turns off walls before evaluation.

Referenced by StatisticsVector< T >::get_walls(), StatisticsVector< T >::set_walls(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::xmaxStat [private]

Referenced by StatisticsVector< T >::get_xmaxStat(), StatisticsVector< T >::set_xmaxStat(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::xminStat [private]

By default the points of evaluation are placed in an grid on the domain [xminStat,xmaxStat]x[yminStat,ymaxStat]x[zminStat,zmaxStat].

By default the domain is set to the domain of the MD problem (indicated by setting the stat-domain values to nan), but can be resized.

Referenced by StatisticsVector< T >::get_xminStat(), StatisticsVector< T >::set_xminStat(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::ymaxStat [private]

Referenced by StatisticsVector< T >::get_ymaxStat(), StatisticsVector< T >::set_ymaxStat(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::yminStat [private]

Referenced by StatisticsVector< T >::get_yminStat(), StatisticsVector< T >::set_yminStat(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::zmaxStat [private]

Referenced by StatisticsVector< T >::get_zmaxStat(), StatisticsVector< T >::set_zmaxStat(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>
Mdouble StatisticsVector< T >::zminStat [private]

Referenced by StatisticsVector< T >::get_zminStat(), StatisticsVector< T >::set_zminStat(), and StatisticsVector< T >::StatisticsVector().

The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 17 2013 13:05:10 for HG-MD by   doxygen 1.7.6.1