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

#include <StatisticsPoint.h>

Inheritance diagram for StatisticsVector< T >:

Public Member Functions
void	constructor ()
	this is the actual constructor, sets up all basic things More...

void	constructor (std::string name)

	StatisticsVector ()
	Basic constructor only calls constructor() More...

	StatisticsVector (std::string name)

	StatisticsVector (StatisticsVector &other)
	Copy constructor. More...

	StatisticsVector (int argc, char *argv[])
	Advanced constructor that accepts arguments from the command line. More...

void	readStatArguments (int argc, char *argv[])

std::string	printStat ()
	Outputs member variable values to a std::string. More...

void	set_statType ()

void	writeOutputFiles ()
	Writes the simulation data onto all the files i.e. .data, .ene, .fstat ... More...

void	print_help ()

void	setNX (int new_)

void	set_hx (Mdouble hx)

void	setNY (int new_)

void	set_hy (Mdouble hy)

void	setNZ (int new_)

void	set_hz (Mdouble hz)

void	setN (int n)

void	set_h (Mdouble h)

int	getNX ()

int	getNY ()

int	getNZ ()

void	setCGTimeMin (Mdouble t)

void	setTimeMaxStat (Mdouble t)

void	setCGTimeAveragingInterval (Mdouble t)

Mdouble	getCGTimeMin ()

Mdouble	getTimeMaxStat ()

Mdouble	getCGTimeAveragingInterval ()

bool	check_current_time_for_statistics ()

void	setCGShape (const char *new_)

void	setCGShape (CG new_)

CG	getCGShape ()

void	setN (int nx_, int ny_, int nz_)

void	getN (int &nx_, int &ny_, int &nz_)

void	setCGWidth (Mdouble w)
	Set CG variables w2 and CG_invvolume. More...

void	setCGWidth2 (Mdouble new_)
	Set CG variables w2 and CG_invvolume. More...

Mdouble	getCGWidth ()

Mdouble	getCGWidthSquared ()

Mdouble	getCutoff ()

Mdouble	getCutoff2 ()

std::string	print_CG ()
	Output coarse graining variables. More...

StatisticsPoint< T >	average (std::vector< StatisticsPoint< T > > &P)
	Output average of statistical variables. More...

virtual void	reset_statistics ()
	Set all statistical variables to zero. More...

void	statistics_from_fstat_and_data ()
	get StatisticsPoint More...

void	statistics_from_p3 ()
	this is a modified version of statistics_from_fstat_and_data. It is used to read p3d files More...

void	jump_p3c ()
	get force statistics from particle collisions More...

void	setDoTimeAverage (bool new_)

bool	getDoTimeAverage ()

void	setStressTypeForFixedParticles (int new_)

int	getStressTypeForFixedParticles ()

void	set_infiniteStressForFixedParticles (bool new_)

void	setMirrorAtDomainBoundary (Mdouble new_)

Mdouble	getMirrorAtDomainBoundary ()

void	setDoVariance (bool new_)

bool	getDoVariance ()

void	setDoGradient (bool new_)

bool	getDoGradient ()

void	setSuperExact (bool new_)

bool	getSuperExact ()

void	setDoIgnoreFixedParticles (bool new_)

bool	getDoIgnoreFixedParticles ()

void	verbose ()

void	setVerbosityLevel (int new_)

int	getVerbosityLevel () const

void	setCGWidthalls (bool new_)

bool	getCGWidthalls ()

void	setDoPeriodicWalls (bool new_)

bool	getDoPeriodicWalls ()

void	setCGWidth_over_rmax (Mdouble new_)

Mdouble	getCGWidth_over_rmax ()

void	setPositions ()
	Set position of StatisticsPoint points and set variables to 0. More...

bool	readNextDataFile (unsigned int format)

void	gather_force_statistics_from_fstat_and_data ()
	get force statistics from particle collisions More...

void	gather_force_statistics_from_p3c (int version)
	get force statistics from particle collisions More...

void	gather_force_statistics_from_p3w (int version, std::vector< int > &index)
	get force statistics from particle collisions More...

void	evaluate_force_statistics (int wp=0)
	get force statistics More...

void	evaluate_wall_force_statistics (Vec3D P, int wp=0)
	get force statistics (i.e. first particle is a wall particle) More...

void	jump_fstat ()

void	initialiseStatistics ()
	Initializes statistics, i.e. setting w2, setting the grid and writing the header lines in the .stat file. More...

void	outputStatistics ()
	Calculates statistics for Particles (i.e. not collisions) More...

void	gatherContactStatistics (unsigned 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) More...

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

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

void	write_statistics ()
	Writes regular statistics. More...

void	write_time_average_statistics ()
	Writes out time averaged statistics. More...

void	evaluate_particle_statistics (std::vector< BaseParticle * >::iterator P, int wp=0)
	Calculates statistics for a single Particle. More...

std::vector< StatisticsPoint< T > >	getCGPoints ()

Mdouble	getXMinStat ()
	Functions to acces and change the domain of statistics. More...

Mdouble	getYMinStat ()

Mdouble	getZMinStat ()

Mdouble	getXMaxStat ()

Mdouble	getYMaxStat ()

Mdouble	getZMaxStat ()

void	setXMinStat (Mdouble xminStat_)

void	setYMinStat (Mdouble yminStat_)

void	setZMinStat (Mdouble zminStat_)

void	setXMaxStat (Mdouble xmaxStat_)

void	setYMaxStat (Mdouble ymaxStat_)

void	setZMaxStat (Mdouble zmaxStat_)

int	getNTimeAverage ()

Mdouble	setInfinitelyLongDistance ()

void	set_Polynomial (std::vector< Mdouble > new_coefficients, unsigned int new_dim)

void	set_Polynomial (Mdouble *new_coefficients, unsigned int num_coeff, unsigned int new_dim)

void	setPolynomialName (const char *new_name)

std::string	getPolynomialName ()

void	setDoDoublePoints (bool new_)

bool	getDoDoublePoints ()

void	setNTimeAverageReset (int new_)

bool	getNTimeAverageReset ()

void	set_rmin (Mdouble new_)

void	set_rmax (Mdouble new_)

void	set_hmax (Mdouble new_)

Mdouble	evaluatePolynomial (Mdouble r)

Mdouble	evaluatePolynomialGradient (Mdouble r)

Mdouble	evaluateIntegral (Mdouble n1, Mdouble n2, Mdouble t)

unsigned int	getStepSize () const

void	setStepSize (unsigned int stepSize)

template<>
void	setNZ (int new_ UNUSED)

template<>
void	setNY (int new_ UNUSED)

template<>
void	setNX (int new_ UNUSED)

template<>
void	setNY (int new_ UNUSED)

template<>
void	setNZ (int new_ UNUSED)

template<>
void	setNX (int new_ UNUSED)

template<>
void	setNZ (int new_ UNUSED)

template<>
void	setNX (int new_ UNUSED)

template<>
void	setNY (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<>
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<>
Mdouble	setInfinitelyLongDistance ()

template<>
Mdouble	setInfinitelyLongDistance ()

template<>
Mdouble	setInfinitelyLongDistance ()

template<>
Mdouble	setInfinitelyLongDistance ()

template<>
Mdouble	setInfinitelyLongDistance ()

template<>
Mdouble	setInfinitelyLongDistance ()

Public Member Functions inherited from DPMBase
void	constructor ()
	A function which initialises the member variables to default values, so that the problem can be solved off the shelf; sets up a basic two dimensional problem which can be solved off the shelf. It is called in the constructor DPMBase(). More...

	DPMBase ()
	Constructor that calls the "void constructor()". More...

	DPMBase (const FilesAndRunNumber &other)
	Copy constructor type-1. More...

	DPMBase (const DPMBase &other)
	Copy constructor type-2. More...

virtual	~DPMBase ()
	virtual destructor More...

void	solve ()
	The work horse of the code. More...

void	checkSettings ()
	Checks if the essentials are set properly to go ahead with solving the problem. More...

void	solve (int argc, char *argv[])
	The solve function is the work horse of the code with the user input. More...

virtual void	setupInitialConditions ()
	This function allows to set the initial conditions for our problem to be solved, by default particle locations are randomly set. Remember particle properties must also be defined here. More...

virtual void	writeXBallsScript () const
	This writes a script which can be used to load the xballs problem to display the data just generated. More...

virtual double	getInfo (const BaseParticle &P) const
	A virtual method that allows the user to overrride and set what is written into the info column in the data file. By default it returns the Species ID number. More...

virtual void	writeRestartFile ()
	Stores all the particle data for current save time step. Calls the write function. More...

int	readRestartFile ()
	Reads all the particle data corresponding to the current saved time step. Which is what the restart file basically stores. The problem description with the latest particle data. More...

int	readRestartFile (std::string fileName)
	Also reads all the particle data corresponding to the current saved time step. More...

virtual void	write (std::ostream &os, bool writeAllParticles=true) const
	Loads all MD data and plots statistics for all timesteps in the .data file. More...

virtual void	read (std::istream &is)
	Reads all particle data into a restart file. More...

virtual void	readOld (std::istream &is)
	Reads all particle data into a restart file; old version. More...

bool	readDataFile (const std::string fileName, unsigned int format=0)
	This allows particle data to be reloaded from data files. More...

bool	readParAndIniFiles (const std::string fileName)
	Allows the user to read par.ini files (useful to read MDCLR files) More...

bool	readNextDataFile (unsigned int format=0)
	Reads the next data file with default format=0. However, one can modify the format based on whether the particle data corresponds to 3D or 2D data. See XBalls/xballs.txt. More...

bool	findNextExistingDataFile (Mdouble tMin, bool verbose=true)
	Useful when fileType is chosen as Multiple Files or Multiple files with padded. More...

bool	readArguments (int argc, char *argv[])
	Can interpret main function input arguments that are passed by the driver codes. More...

virtual bool	readNextArgument (int &i, int argc, char *argv[])

virtual bool	checkParticleForInteraction (const BaseParticle &P)
	Checks if the particle having any interaction with walls or other particles. More...

Mdouble	getTime () const
	Access function for the time. More...

unsigned int	getNtimeSteps () const
	Returns the current counter of time steps. More...

void	setTime (Mdouble time)
	Access function for the time. More...

void	setTimeMax (Mdouble newTMax)
	Allows the upper time limit to be changed. More...

Mdouble	getTimeMax () const
	Allows the user to access the total simulation time during the simulation. Cannot change it though. More...

void	setDoCGAlways (bool newDoCGFlag)

void	setRotation (bool newRotFlag)
	Allows to set the flag for enabling or disabling particle rotation in the simulations. More...

bool	getRotation () const
	Returns a flag indicating if particle rotation is enabled or disabled. More...

bool	getDoCGAlways () const

Mdouble	getXMin () const
	If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin,. More...

Mdouble	getXMax () const
	If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax,. More...

Mdouble	getYMin () const
	If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin,. More...

Mdouble	getYMax () const
	If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax,. More...

Mdouble	getZMin () const
	If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin,. More...

Mdouble	getZMax () const
	If the length of the problem domain in z-direction is ZMax - ZMin, then getZMax() returns ZMax,. More...

void	setXMin (Mdouble newXMin)
	If the length of the problem domain in x-direction is XMax - XMin, this method sets XMin. More...

void	setYMin (Mdouble newYMin)
	If the length of the problem domain in y-direction is YMax - YMin, this method sets YMin. More...

void	setZMin (Mdouble newZMin)
	If the length of the problem domain in z-direction is ZMax - ZMin, this method sets ZMin. More...

void	setXMax (Mdouble newXMax)
	If the length of the problem domain in x-direction is XMax - XMin, this method sets XMax. More...

void	setYMax (Mdouble newYMax)
	If the length of the problem domain in y-direction is YMax - YMin, this method sets YMax. More...

void	setZMax (Mdouble newZMax)
	If the length of the problem domain in z-direction is XMax - XMin, this method sets ZMax. More...

void	setTimeStep (Mdouble newDt)
	Allows the time step dt to be changed. More...

Mdouble	getTimeStep () const
	Allows the time step dt to be accessed. More...

void	setXBallsColourMode (int newCMode)
	Set the xball output mode. More...

int	getXBallsColourMode () const
	Get the xball colour mode (CMode) More...

void	setXBallsVectorScale (double newVScale)
	Set the scale of vectors in xballs. More...

double	getXBallsVectorScale () const
	Returns the scale of vectors used in xballs. More...

void	setXBallsAdditionalArguments (std::string newXBArgs)
	Set the additional arguments for xballs. More...

std::string	getXBallsAdditionalArguments () const

void	setXBallsScale (Mdouble newScale)
	Sets the scale of the view (either normal, zoom in or zoom out) to display in xballs. The default is fit to screen. More...

double	getXBallsScale () const
	Returns the scale of the view in xballs. More...

void	setGravity (Vec3D newGravity)
	Allows to modify the gravity vector. More...

Vec3D	getGravity () const
	Returns the gravity vector. More...

void	setDimension (unsigned int newDim)
	Sets the system and particle dimension. More...

void	setSystemDimensions (unsigned int newDim)
	Allows for the dimension of the simulation to be changed. More...

unsigned int	getSystemDimensions () const
	Returns the dimension of the simulation. Note there is also a particle dimension. More...

void	setParticleDimensions (unsigned int particleDimensions)
	Allows the dimension of the particle (f.e. for mass) to be changed. e.g. discs or spheres. More...

unsigned int	getParticleDimensions () const
	Returns the particle dimensions. More...

std::string	getRestartVersion () const
	This is to take into account for different Mercury versions. Returns the version of the restart file. More...

void	setRestartVersion (std::string newRV)
	Sets restart_version. More...

bool	getRestarted () const
	Returns the flag denoting if the simulation was restarted or not. More...

void	setRestarted (bool newRestartedFlag)
	Allows to set the flag stating if the simulation is to be restarted or not. More...

bool	getAppend () const
	Returns the flag denoting if the append option is on or off. More...

void	setAppend (bool newAppendFlag)
	Allows to set the append option. More...

Mdouble	getElasticEnergy () const
	Returns the global elastic energy within the system. More...

Mdouble	getKineticEnergy () const
	Returns the global kinetic energy stored in the system. More...

bool	areInContact (const BaseParticle pI, const BaseParticle pJ) const
	Checks if two particle are in contact or is there any positive overlap. More...

virtual void	hGridInsertParticle (BaseParticle *obj UNUSED)
	no implementation but can be overidden in its derived classes. More...

virtual void	hGridUpdateParticle (BaseParticle *obj UNUSED)
	no implementation but can be overidden in its derived classes. More...

virtual void	hGridRemoveParticle (BaseParticle *obj UNUSED)
	no implementation but can be overidden in its derived classes. More...

virtual void	hGridUpdateMove (BaseParticle *, Mdouble)
	no implementation but can be overidden in its derived classes. More...

Public Member Functions inherited from FilesAndRunNumber
	FilesAndRunNumber ()
	Constructor. More...

	FilesAndRunNumber (const FilesAndRunNumber &other)
	Copy constructor. More...

virtual	~FilesAndRunNumber ()
	Constructor. More...

void	constructor ()
	a function called by the FilesAndRunNumber() (constructor) More...

void	incrementRunNumberInFile ()
	Increment the run Number (counter value) stored in the file_counter (COUNTER_DONOTDEL) by 1 and store the new value in the counter file. More...

int	readRunNumberFromFile ()
	Read the run number or the counter from the counter file (COUNTER_DONOTDEL) More...

void	autoNumber ()
	The autoNumber() function is the trigger. It calls three functions. setRunNumber(), readRunNumberFromFile() and incrementRunNumberInFile(). More...

std::vector< int >	get2DParametersFromRunNumber (int size_x, int size_y)
	This turns a counter into two indices which is an amazing feature for doing two dimensional parameter studies. The indices run from 1:size_x and 1:size_y, while the study number starts at 0 ( initially the counter=1 in COUNTER_DONOTDEL) More...

int	launchNewRun (const char *name, bool quick=false)
	This launches a code from within this code. Please pass the name of the code to run. More...

void	setRunNumber (int runNumber)
	This sets the counter/Run number, overriding the defaults. More...

int	getRunNumber () const
	This returns the current value of the counter (runNumber_) More...

void	read (std::istream &is)
	Accepts an input stream std::istream. More...

void	write (std::ostream &os) const
	Accepts an output stream read function, which accepts an input stream std::ostream. More...

Public Member Functions inherited from Files
	Files ()
	A constructor. More...

virtual	~Files ()
	A destructor, watch out its a virtual destructor. More...

	Files (const Files &other)
	Copy constructor. More...

MERCURY_DEPRECATED File &	getDataFile ()
	The non const version. Allows one to edit the File::dataFile. More...

MERCURY_DEPRECATED File &	getEneFile ()
	The non const version. Allows to edit the File::eneFile. More...

MERCURY_DEPRECATED File &	getFStatFile ()
	The non const version. Allows to edit the File::fStatFile. More...

MERCURY_DEPRECATED File &	getRestartFile ()
	The non const version. Allows to edit the File::restartFile. More...

MERCURY_DEPRECATED File &	getStatFile ()
	The non const version. Allows to edit the File::statFile. More...

MERCURY_DEPRECATED const File &	getDataFile () const
	The const version. Does not allow for any editing of the File::dataFile. More...

MERCURY_DEPRECATED const File &	getEneFile () const
	The const version. Does not allow for any editing of the File::eneFile. More...

MERCURY_DEPRECATED const File &	getFStatFile () const
	The const version. Does not allow for any editing of the File::fStatFile. More...

MERCURY_DEPRECATED const File &	getRestartFile () const
	The const version. Does not allow for any editing of the File::restartFile. More...

MERCURY_DEPRECATED const File &	getStatFile () const
	The const version. Does not allow for any editing of the File::statFile. More...

const std::string &	getName () const
	Returns the name of the file. Does not allow to change it though. More...

void	setName (const std::string &name)
	Allows to set the name of all the files (ene, data, fstat, restart, stat) More...

void	setName (const char *name)
	Calls setName(std::string) More...

void	setSaveCount (unsigned int saveCount)
	Sets File::saveCount_ for all files (ene, data, fstat, restart, stat) More...

void	setFileType (FileType fileType)
	Sets File::fileType_ for all files (ene, data, fstat, restart, stat) More...

void	setOpenMode (std::fstream::openmode openMode)
	Sets File::openMode_ for all files (ene, data, fstat, restart, stat) More...

void	resetFileCounter ()
	Resets the file counter for each file i.e. for ene, data, fstat, restart, stat) More...

void	read (std::istream &is)
	Extracts data from the input stream (which is basically a file you want to read from) into name_, restartFile .... More...

void	write (std::ostream &os) const
	Writes data into a file from the member variables name_, restartFile, dataFile etc. More...

void	openFiles ()
	Opens all the files (ene, data, fstat, restart, stat) for reading and writing purposes. More...

void	closeFiles ()
	Closes all files (ene, data, fstat, restart, stat) that were opened to read or write. More...

void	setNextSavedTimeStep (unsigned int nextSavedTimeStep)
	Sets the next time step for all the files (ene, data, fstat, restart, stat) at which the data is to be written or saved. More...

Protected Member Functions
bool	satisfiesInclusionCriteria (BaseParticle *P)

bool	loadVelocityProfile (const char *filename)

bool	loadPositions (const char *filename)

Vec3D	getVelocityProfile (Vec3D Position)

bool	read_next_from_p3p_file ()

void	auto_setdim ()

Protected Member Functions inherited from DPMBase
virtual void	computeAllForces ()
	Computes all the forces acting on the particles by using the setTorque and setForce methods. See BaseInteractible.cc. More...

virtual void	computeInternalForces (BaseParticle *i)
	Computes the forces between particles (internal in the sense that the sum over all these forces is zero i.e. fully modelled forces) More...

virtual void	computeInternalForces (BaseParticle P1, BaseParticle P2)
	Computes the forces between particles (internal in the sense that the sum over all these forces is zero i.e. fully modelled forces) More...

virtual void	computeExternalForces (BaseParticle *PI)
	Computes the external forces acting on particles (e.g. gravitational) More...

virtual void	computeForcesDueToWalls (BaseParticle *PI)
	Computes the forces on the particles due to the walls (normals are outward normals) More...

virtual void	actionsOnRestart ()
	A virtual function where the users can add extra code which is executed only when the code is restarted. More...

virtual void	actionsBeforeTimeLoop ()
	A virtual function. Allows one to carry out any operations before the start of the time loop. More...

virtual void	hGridActionsBeforeTimeLoop ()
	A virtual function that allows one to carry out hGrid operations before the start of the time loop. More...

virtual void	hGridActionsBeforeTimeStep ()
	A virtual function that allows one to set or execute hGrid parameters or operations before every simulation time step. More...

virtual bool	getHGridUpdateEachTimeStep () const

virtual void	actionsBeforeTimeStep ()
	A virtual function which allows to define operations to be executed before the new time step. More...

virtual void	actionsAfterSolve ()
	A virtual function which allows to define operations to be executed after the solve(). More...

virtual void	actionsAfterTimeStep ()
	A virtual function which allows to define operations to be executed after time step. More...

virtual void	outputXBallsData (std::ostream &os) const
	This function writes the location of the walls and particles in a format the XBalls program can read. See XBalls/xballs.txt. However, MercuryDPM supports a much better viewer now called Paraview. See the tutorials section in the documentation. More...

virtual void	outputXBallsDataParticle (const unsigned int i, const unsigned int format, std::ostream &os) const
	This function writes out the particle locations into an output stream in a format the XBalls program can read. More...

virtual void	writeEneHeader (std::ostream &os) const
	Writes a header with a certain format for ENE file. More...

virtual void	writeFstatHeader (std::ostream &os) const
	Writes a header with a certain format for FStat file. More...

virtual void	writeEneTimestep (std::ostream &os) const
	This function enables one to write the global energy available in the system after each time step. The default is to compute the rotational and translational kinetic energy, potential energy and the centre of mass. More...

void	gatherContactStatistics ()

virtual void	integrateBeforeForceComputation ()
	This is were the integration is done, at the moment it is velocity Verlet integration and is done before the forces are computed. See http://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet. More...

virtual void	integrateAfterForceComputation ()
	Integration is done after force computations. We apply the Velocity verlet scheme. See http://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet. More...

virtual void	checkInteractionWithBoundaries ()
	There are a range of boundaries one could implement depending on ones' problem. This methods checks for interactions between particles and such range of boundaries. See BaseBoundary.h and all the boundaries in the Boundaries folder. More...

virtual void	hGridActionsBeforeIntegration ()
	no implementation but can be overidden in its derived classes. More...

virtual void	hGridActionsAfterIntegration ()
	no implementation but can be overidden in its derived classes. More...

virtual void	broadPhase (BaseParticle *i)
	By broad one means to screen and determine an approximate number of particles that a given particle can be in contact with. Hence the word "Broad phase" of contact detection. More...

void	setFixedParticles (unsigned int n)

void	initialiseTangentialSprings ()

virtual void	printTime () const
	Displays the current simulation time onto your screen for example. More...

virtual bool	continueSolve () const

void	outputInteractionDetails () const
	Displays the interaction details corresponding to the pointer objects in the interaction handler. More...

bool	isTimeEqualTo (Mdouble time) const
	Checks if the input variable "time" is the current time in the simulation. More...

void	removeDuplicatePeriodicParticles ()
	Removes periodic duplicate Particles. More...

void	checkAndDuplicatePeriodicParticles ()
	In case of periodic boundaries, the below methods checks and adds particles when necessary into the particle handler. See DPMBase.cc and PeriodicBoundary.cc for more details. More...

Protected Attributes
StatType	statType
	Possible values X,Y,Z,XY,XZ,YZ,XYZ,RAZ,RA,RZ,AZ,R,A are used to determine if the statistics are averaged; f.e. StatType X is averaged over y and z. More...

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]. More...

int	ny
	see nx More...

int	nz
	see nx More...

Mdouble	xminStat
	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. More...

Mdouble	xmaxStat

Mdouble	yminStat

Mdouble	ymaxStat

Mdouble	zminStat

Mdouble	zmaxStat

int	nxMirrored
	extension of grid size from mirrored points More...

int	nyMirrored

int	nzMirrored

std::vector< StatisticsPoint< T > >	Points
	A vector that stores the values of the statistical variables at a given position. More...

std::vector< StatisticsPoint< T > >	dx
	A vector that stores the gradient in x of all statistical variables at a given position. More...

std::vector< StatisticsPoint< T > >	dy
	A vector that stores the gradient in y of all statistical variables at a given position. More...

std::vector< StatisticsPoint< T > >	dz
	A vector that stores the gradient in z of all statistical variables at a given position. More...

std::vector< StatisticsPoint< T > >	timeAverage
	A vector used to sum up all statistical values in Points for time-averaging. More...

std::vector< StatisticsPoint< T > >	timeVariance
	a vector used to sum up the variance in time of all statistical values More...

std::vector< StatisticsPoint< T > >	dxTimeAverage
	a vector used to sum up all statistical gradients in dx for time-averaging More...

std::vector< StatisticsPoint< T > >	dyTimeAverage
	a vector used to sum up all statistical gradients in dy for time-averaging More...

std::vector< StatisticsPoint< T > >	dzTimeAverage
	a vector used to sum up all statistical gradients in dz for time-averaging More...

bool	doTimeAverage
	Determines if output is averaged over time. More...

int	nTimeAverageReset
	Determines after how many timesteps the time average is reset. More...

bool	doVariance
	Determines if variance is outputted. More...

bool	doGradient
	Determines if gradient is outputted. More...

int	nTimeAverage
	A counter needed to average over time. More...

CG	CG_type
	coarse graining type (Gaussian, Heaviside, Polynomial) More...

NORMALIZED_POLYNOMIAL< T >	CGPolynomial
	Stores the Polynomial, if the cg function is an axisymmetric function Polynomial in r. More...

Mdouble	w2
	coarse graining width squared; for HeavisideSphere and Gaussian More...

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

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 More...

Mdouble	tminStat
	Statistical output will only be created if t>tminStat. More...

Mdouble	tmaxStat
	Statistical output will only be created if t<tmaxStat. More...

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

Mdouble	indSpecies
	defines the species for which statistics are extracted (-1 for all species) More...

Mdouble	rmin
	defines the minimum radius of the particles for which statistics are extracted More...

Mdouble	rmax
	defines the maximum radius of the particles for which statistics are extracted More...

Mdouble	hmax
	defines the maximum height of the particles for which statistics are extracted More...

bool	walls
	Turns off walls before evaluation. More...

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

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

int	StressTypeForFixedParticles
	Stress type for fixed particles. More...

int	verbosity
	Determines how much is outputted to the terminal. More...

int	format
	format of the data input file More...

Mdouble	mirrorAtDomainBoundary
	0: Statistics near the wall are equal to statistics away from the wall; 1: Statistics are mirrored at the domain boundaries; up to a maximum depth specified by this number More...

bool	isMDCLR

bool	superexact
	If true, cutoff radius for Gaussian is set to 5w (from 3w) More...

bool	doDoublePoints

Vec3D	P1
	Position of first contact point. More...

Vec3D	P2
	Position of second contact point. More...

Vec3D	P1_P2_normal
	Direction of contact. More...

Mdouble	P1_P2_distance
	Length of contact line. More...

Matrix3D	P1_P2_NormalStress
	Contact stress from normal forces along the line of contact. More...

Matrix3D	P1_P2_ContactCoupleStress

Vec3D	P1_P2_Contact

Matrix3D	P1_P2_TangentialStress
	Contact stress from tangential forces along the line of contact. More...

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

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

MatrixSymmetric3D	P1_P2_Fabric
	Fabric. More...

Vec3D	P1_P2_CollisionalHeatFlux
	not yet working More...

Mdouble	P1_P2_Dissipation
	not yet working More...

Mdouble	P1_P2_Potential
	not yet working More...

std::vector< Vec3D >	VelocityProfile

Vec3D	VelocityProfile_Min

Vec3D	VelocityProfile_D

std::fstream	p3p_file

std::fstream	p3c_file

std::fstream	p3w_file

unsigned int	stepSize_

std::vector< Vec3D >	positions_

Additional Inherited Members
Public Attributes inherited from DPMBase
SpeciesHandler	speciesHandler
	A handler to that stores the species type i.e. elastic, linear visco-elastic... et cetera. More...

RNG	random
	This is a random generator, often used for setting up the initial conditions etc... More...

ParticleHandler	particleHandler
	An object of the class ParticleHandler, contains the pointers to all the particles created. More...

WallHandler	wallHandler
	An object of the class WallHandler. Contains pointers to all the walls created. More...

BoundaryHandler	boundaryHandler
	An object of the class BoundaryHandler which concerns insertion and deletion of particles into or from regions. More...

InteractionHandler	interactionHandler
	An object of the class InteractionHandler. More...

Public Attributes inherited from Files
File	dataFile
	An instance of class File to handle in- and output into a .data file. More...

File	fStatFile
	An instance of class File to handle in- and output into a .fstat file. More...

File	eneFile
	An instance of class File to handle in- and output into a .ene file. More...

File	restartFile
	An instance of class File to handle in- and output into a .restart file. More...

File	statFile
	An instance of class File to handle in- and output into a .stat file. More...

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).

Definition at line 36 of file StatisticsPoint.h.

Constructor & Destructor Documentation

template<StatType T>

StatisticsVector< T >::StatisticsVector ( )

inline

Basic constructor only calls constructor()

Definition at line 75 of file StatisticsVector.h.

     {
         constructor();
     }

template<StatType T>

StatisticsVector< T >::StatisticsVector ( std::string name )

inline

Definition at line 83 of file StatisticsVector.h.

     {
         constructor(name);
     }

template<StatType T>

StatisticsVector< T >::StatisticsVector ( StatisticsVector< T > & other )

Copy constructor.

Definition at line 180 of file StatisticsVector.hcc.

         : DPMBase()
 {
     //assumes that we want the statistical method copied, but resets the time averaging 
     constructor();
     //~ StatType statType; //this has to be constant, right?
     nx = other.nx;
     ny = other.ny;
     nz = other.nz;
     xminStat = other.xminStat;
     yminStat = other.yminStat;
     zminStat = other.zminStat;
     xmaxStat = other.xmaxStat;
     ymaxStat = other.ymaxStat;
     zmaxStat = other.zmaxStat;
     nxMirrored = other.nxMirrored;
     nyMirrored = other.nyMirrored;
     nzMirrored = other.nzMirrored;
     Points = other.Points;
     dx = other.dx;
     dy = other.dy;
     dz = other.dz;
     
     timeAverage = other.timeAverage;
     timeVariance = other.timeVariance;
     dxTimeAverage = other.dxTimeAverage;
     dyTimeAverage = other.dyTimeAverage;
     dzTimeAverage = other.dzTimeAverage;
     doTimeAverage = other.doTimeAverage;
     nTimeAverage = 0;
     for (unsigned int i = 0; i < timeAverage.size(); i++)
         timeAverage[i].set_zero();
     for (unsigned int i = 0; i < timeVariance.size(); i++)
         timeVariance[i].set_zero();
     for (unsigned int i = 0; i < dxTimeAverage.size(); i++)
         dxTimeAverage[i].set_zero();
     for (unsigned int i = 0; i < dyTimeAverage.size(); i++)
         dyTimeAverage[i].set_zero();
     for (unsigned int i = 0; i < dzTimeAverage.size(); i++)
         dzTimeAverage[i].set_zero();
     doVariance = other.doVariance;
     doGradient = other.doGradient;
     CG_type = other.CG_type;
     CGPolynomial = other.CGPolynomial;
     w2 = other.w2;
     cutoff = other.cutoff;
     cutoff2 = other.cutoff2;
     w_over_rmax = other.w_over_rmax;
     tminStat = other.tminStat;
     tmaxStat = other.tmaxStat;
     tintStat = other.tintStat;
     rmin = other.rmin;
     rmax = other.rmax;
     hmax = other.hmax;
     indSpecies = other.indSpecies;
     walls = other.walls;
     periodicWalls = other.periodicWalls;
     ignoreFixedParticles = other.ignoreFixedParticles;
     StressTypeForFixedParticles = other.StressTypeForFixedParticles;
     verbosity = other.verbosity;
     format = other.format;
     mirrorAtDomainBoundary = other.mirrorAtDomainBoundary;
     isMDCLR = other.isMDCLR;
     superexact = other.superexact;
 }

template<StatType T>

StatisticsVector< T >::StatisticsVector	(	int	argc,
		char *	argv[]
	)

Advanced constructor that accepts arguments from the command line.

Definition at line 279 of file StatisticsVector.hcc.

 {
     //check if filename is given
     if (argc < 2)
     {
         std::cerr << "fstatistics.exe filename [-options]" << std::endl
                 << "Type -help for more information" << std::endl;
         exit(-1);
     }
     
     //print help if required
     if (!strcmp(argv[1], "-help"))
     {
         print_help();
     }
     
     //create StatisticsVector
     std::string filename(argv[1]);
     constructor(filename);
 
     //check for options (standardly '-option argument')
     readStatArguments(argc, argv);
 }

Member Function Documentation

template<StatType T>

void StatisticsVector< T >::auto_setdim ( )

protected

Definition at line 1243 of file StatisticsVector.hcc.

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

template<StatType T>

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

Output average of statistical variables.

Output names of statistical variables.

Definition at line 872 of file StatisticsVector.hcc.

 {
     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

Definition at line 269 of file StatisticsVector.h.

     {
         return (getTime() >= getCGTimeMin() && getTime() <= getTimeMaxStat() + getTimeStep());
     }

template<StatType T>

void StatisticsVector< T >::constructor ( )

this is the actual constructor, sets up all basic things

Definition at line 115 of file StatisticsVector.hcc.

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

Referenced by StatisticsVector< O >::StatisticsVector(), and 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;
     setCGTimeAveragingInterval(0);
     setCGTimeMin(-1e20);
     setTimeMaxStat(NAN);
     rmin = 0;
     hmax = 1e20;
     rmax = 1e20;
     indSpecies = -1; //all species
             
     //set default options - Not used it physially don't make sense i.e if ?max<?min
     xminStat = NAN;
     xmaxStat = NAN;
     yminStat = NAN;
     ymaxStat = NAN;
     zminStat = NAN;
     zmaxStat = NAN;
     
     // set default options
     ignoreFixedParticles = true; //this should be true if you want statistics only of the flow
     StressTypeForFixedParticles = 1;
     verbosity = 1;
     walls = true;
     periodicWalls = true;
     format = 0;
     mirrorAtDomainBoundary = false;
     setSuperExact(false);
     VelocityProfile.resize(0); //equivalent to do not use velocityprofile
     setStepSize(1);
     
     // time averaging
     setDoTimeAverage(true);
     nTimeAverage = 0;
     nTimeAverageReset = -1;
     //calculate variance
     setDoVariance(false);
     //calculate gradient
     setDoGradient(false);
     doDoublePoints = false;
     
     // additional stuff
     statFile.setFileType(FileType::ONE_FILE);
 }

template<StatType T>

void StatisticsVector< T >::constructor ( std::string name )

Todo:: TWnow do we really need a default species?

Definition at line 247 of file StatisticsVector.hcc.

 {
     //call default constructor
     constructor();
     
     // set name
     setName(name.c_str());
 
     if (!strcmp(getName().c_str(), "c3d"))
     {
         //write here what to do with Stefan's data 
         std::cout << "MDCLR data" << std::endl;
         exit(-1);
     }
     else
     {
         //write here what to do with Mercury data
         isMDCLR = false;
         if (readRestartFile()==0)
         {
             std::cout << "no restart file given: using default density 1" << std::endl;
             speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
 
         }
         setAppend(false);
         //~ setTimeMaxStat(t+dt); //note:doesn't work if restart data is not the latest
     }
 }

template<StatType T>

void StatisticsVector< T >::evaluate_force_statistics ( int wp = 0 )

get force statistics

todo{What are P1 and P2 over here?}

todo{DK: This only works for PeriodicBoundary and not for AngledPeriodicBoundaries}

Todo:: {TW: Couple stress calculation works only for Gaussians}

Definition at line 2207 of file StatisticsVector.hcc.

References Matrix3D::cross(), AngledPeriodicBoundary::distance(), PeriodicBoundary::getDistance(), AngledPeriodicBoundary::shiftPositions(), PeriodicBoundary::shiftPositions(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

 {
     if (periodicWalls)
     {
         for (unsigned int k = wp; k < boundaryHandler.getNumberOfObjects(); k++)
         {
             PeriodicBoundary *b0 = dynamic_cast<PeriodicBoundary *>(boundaryHandler.getObject(k));
             if (b0)
             {
                 b0->getDistance(P1);
                 b0->shiftPositions(P1, P2);
                 evaluate_force_statistics(k + 1);
                 b0->shiftPositions(P1, P2);
             }
             AngledPeriodicBoundary* b1 = dynamic_cast<AngledPeriodicBoundary*>(boundaryHandler.getObject(k));
             if (b1 != nullptr)
             {
                 b1->distance(P1);
                 b1->shiftPositions(P1, P2);
                 evaluate_force_statistics(k + 1);
                 b1->shiftPositions(P1, P2);
             }
         }
     }
 
     //evaluate fields that only depend on CParticle parameters
     Mdouble psi; //Course graining integral
     Vec3D rpsi = Vec3D(0, 0, 0);
     static unsigned int counter = 0;
     for (unsigned int i = 0; i < Points.size(); i++)
     {
         psi = Points[i].CG_integral(P1, P2, P1_P2_normal, P1_P2_distance, rpsi);
         if (psi != 0.0)
         {
             counter++;
             if (!std::isfinite(psi))
             {
                 std::cerr << "error: psi =" << psi << " is infinite for "
                     <<   "P1=" << P1
                     << ", P2=" << P2
                     << ", P1_P2_normal=" << P1_P2_normal
                     << ", P1_P2_distance=" << P1_P2_distance
                     << ", t=" << getTime() << std::endl;
                 psi = 0;
             }
             Points[i].ContactCoupleStress += (Matrix3D::cross(P1_P2_Contact * psi, P1_P2_ContactCoupleStress) - Matrix3D::cross(rpsi, P1_P2_ContactCoupleStress)) * (-0.5);
             Points[i].NormalStress += P1_P2_NormalStress * psi;
             Points[i].TangentialStress += P1_P2_TangentialStress * psi;
             Points[i].Fabric += P1_P2_Fabric * psi;
             Points[i].CollisionalHeatFlux += P1_P2_CollisionalHeatFlux * psi;
             Points[i].Dissipation += P1_P2_Dissipation * psi;
             Points[i].Potential += P1_P2_Potential * psi;
             if (getDoGradient())
             {
                 Vec3D d_psi = Points[i].CG_integral_gradient(P1, P2, P1_P2_normal, P1_P2_distance);
                 
                 dx[i].NormalStress += P1_P2_NormalStress * d_psi.X;
                 dx[i].TangentialStress += P1_P2_TangentialStress * d_psi.X;
                 dx[i].Fabric += P1_P2_Fabric * d_psi.X;
                 dx[i].CollisionalHeatFlux += P1_P2_CollisionalHeatFlux * d_psi.X;
                 dx[i].Dissipation += P1_P2_Dissipation * d_psi.X;
                 dx[i].Potential += P1_P2_Potential * d_psi.X;
                 
                 dy[i].NormalStress += P1_P2_NormalStress * d_psi.Y;
                 dy[i].TangentialStress += P1_P2_TangentialStress * d_psi.Y;
                 dy[i].Fabric += P1_P2_Fabric * d_psi.Y;
                 dy[i].CollisionalHeatFlux += P1_P2_CollisionalHeatFlux * d_psi.Y;
                 dy[i].Dissipation += P1_P2_Dissipation * d_psi.Y;
                 dy[i].Potential += P1_P2_Potential * d_psi.Y;
                 
                 dz[i].NormalStress += P1_P2_NormalStress * d_psi.Z;
                 dz[i].TangentialStress += P1_P2_TangentialStress * d_psi.Z;
                 dz[i].Fabric += P1_P2_Fabric * d_psi.Z;
                 dz[i].CollisionalHeatFlux += P1_P2_CollisionalHeatFlux * d_psi.Z;
                 dz[i].Dissipation += P1_P2_Dissipation * d_psi.Z;
                 dz[i].Potential += P1_P2_Potential * d_psi.Z;
             } //end if getDoGradient
         } //end if phi
     } // end forall Points
       //std::cout << "NCeval" << counter << std::endl;
 }

template<StatType T>

void StatisticsVector< T >::evaluate_particle_statistics	(	std::vector< BaseParticle * >::iterator	P,
		int	wp = `0`
	)

Calculates statistics for a single Particle.

Todo:: {I currently disabled displacement calculation as it takes forever}

begin: Linear displacement, $1/(2 \rho^2) \sum_{pq} m_p m_q \phi_q *(v_{pqa} \partial_b \phi_p + v_{pqb} \partial_a \phi_p)$

Definition at line 2365 of file StatisticsVector.hcc.

References Vec3D::cross(), AngledPeriodicBoundary::distance(), Matrix3D::dyadic(), PeriodicBoundary::getDistance(), DPMBase::getTime(), MatrixSymmetric3D::selfDyadic(), AngledPeriodicBoundary::shiftPosition(), PeriodicBoundary::shiftPosition(), mathsFunc::square(), MatrixSymmetric3D::symmetrisedDyadic(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

 {
     bool doDisplacement = false;
 
     if (periodicWalls)
     {
         for (unsigned int k = wp; k < boundaryHandler.getNumberOfObjects(); k++)
         {
             PeriodicBoundary* b = dynamic_cast<PeriodicBoundary*>(boundaryHandler.getObject(k));
             if (b != nullptr)
             {
                 Mdouble dist = b->getDistance(P1);
                 if (mathsFunc::square(dist - (*P)->getInteractionRadius()) < getCutoff2())
                 {
                     b->shiftPosition(*P);
                     evaluate_particle_statistics(P, k + 1);
                     b->shiftPosition(*P);
                 }
             }
             AngledPeriodicBoundary* c = dynamic_cast<AngledPeriodicBoundary*>(boundaryHandler.getObject(k));
             if (c != nullptr)
             {
                 Mdouble dist = c->distance(P1);
                 if (mathsFunc::square(dist - (*P)->getInteractionRadius()) < getCutoff2())
                 {
                     c->shiftPosition(*P);
                     evaluate_particle_statistics(P, k + 1);
                     c->shiftPosition(*P);
                 }
             }
         }
     }
     
     //evaluate fields that only depend on CParticle parameters
     Mdouble phi; //Course graining function
     Vec3D V = Vec3D(0, 0, 0);
     if (VelocityProfile.size())
         V = getVelocityProfile((*P)->getPosition());
     
     for (unsigned int i = 0; i < Points.size(); i++)
     {
         phi = Points[i].CG_function((*P)->getPosition());
         if (phi != 0.0)
         {
             if (!std::isfinite(phi))
             {
                 std::cout
                 << "t =" << DPMBase::getTime()
                         << "error: phi =" << phi
                         << " is infinite for P=" << (*P)->getPosition()
                         << ", Point=" << Points[i].getPosition()
                         << std::endl;
                 phi = 0;
             }
             Points[i].Nu += (*P)->getVolume() * phi;
             Points[i].Density += (*P)->getMass() * phi;
             if (VelocityProfile.size())
             {
                 Points[i].Momentum += ((*P)->getVelocity() - V) * ((*P)->getMass() * phi);
                 Points[i].MomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getVelocity() - V) * ((*P)->getMass() * phi);
             }
             else
             {
                 Points[i].Momentum += (*P)->getVelocity() * ((*P)->getMass() * phi);
                 Points[i].MomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getVelocity()) * ((*P)->getMass() * phi);
             }
             Points[i].DisplacementMomentum += (*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()) * ((*P)->getMass() * phi);
             Points[i].DisplacementMomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount())) * ((*P)->getMass() * phi);
             Points[i].EnergyFlux += (*P)->getVelocity() * ((*P)->getMass() * (*P)->getVelocity().getLengthSquared() / 2 * phi);
             
             Vec3D LocalAngularMomentum = Vec3D::cross(Points[0].clearAveragedDirections((*P)->getPosition() - Points[i].getPosition()), (*P)->getVelocity()) * (*P)->getMass()
                     + (*P)->getAngularVelocity() * (*P)->getInertia();
             Points[i].LocalAngularMomentum += phi * LocalAngularMomentum;
             Points[i].LocalAngularMomentumFlux += Matrix3D::dyadic(LocalAngularMomentum, (*P)->getVelocity() * phi);
             
             //Note: Displacement is only computed directly if gradients are cmputed
             if (getDoGradient())
             {
                 
                 Vec3D d_phi = Points[i].CG_gradient((*P)->getPosition(), phi);
                 
                 if (doDisplacement)
                 {
                     for (std::vector<BaseParticle*>::iterator Q = particleHandler.begin(); Q != particleHandler.end(); ++Q)
                     {
                         double phiQ = Points[i].CG_function((*Q)->getPosition()); //Course graining function
                         Points[i].Displacement += MatrixSymmetric3D::symmetrisedDyadic((*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()) - (*Q)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()), d_phi) * ((*P)->getMass() * (*Q)->getMass() * phiQ);
                     }
                     //end:Displacement
                 }
                 
                 dx[i].Nu += (*P)->getVolume() * d_phi.X;
                 dx[i].Density += (*P)->getMass() * d_phi.X;
                 dx[i].Momentum += (*P)->getVelocity() * ((*P)->getMass() * d_phi.X);
                 dx[i].DisplacementMomentum += (*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()) * ((*P)->getMass() * d_phi.X);
                 dx[i].MomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getVelocity()) * ((*P)->getMass() * d_phi.X);
                 dx[i].DisplacementMomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount())) * ((*P)->getMass() * d_phi.X);
                 dx[i].EnergyFlux += (*P)->getVelocity() * ((*P)->getMass() * (*P)->getVelocity().getLengthSquared() / 2 * d_phi.X);
                 
                 dy[i].Nu += (*P)->getVolume() * d_phi.Y;
                 dy[i].Density += (*P)->getMass() * d_phi.Y;
                 dy[i].Momentum += (*P)->getVelocity() * ((*P)->getMass() * d_phi.Y);
                 dy[i].DisplacementMomentum += (*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()) * ((*P)->getMass() * d_phi.Y);
                 dy[i].MomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getVelocity()) * ((*P)->getMass() * d_phi.Y);
                 dy[i].DisplacementMomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount())) * ((*P)->getMass() * d_phi.Y);
                 dy[i].EnergyFlux += (*P)->getVelocity() * ((*P)->getMass() * (*P)->getVelocity().getLengthSquared() / 2 * d_phi.Y);
                 
                 dz[i].Nu += (*P)->getVolume() * d_phi.Z;
                 dz[i].Density += (*P)->getMass() * d_phi.Z;
                 dz[i].Momentum += (*P)->getVelocity() * ((*P)->getMass() * d_phi.Z);
                 dz[i].DisplacementMomentum += (*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount()) * ((*P)->getMass() * d_phi.Z);
                 dz[i].MomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getVelocity()) * ((*P)->getMass() * d_phi.Z);
                 dz[i].DisplacementMomentumFlux += MatrixSymmetric3D::selfDyadic((*P)->getDisplacement2(getXMin(), getXMax(), getYMin(), getYMax(), getZMin(), getZMax(), getTimeStep() * dataFile.getSaveCount())) * ((*P)->getMass() * d_phi.Z);
                 dz[i].EnergyFlux += (*P)->getVelocity() * ((*P)->getMass() * (*P)->getVelocity().getLengthSquared() / 2 * d_phi.Z);
             } //end if getDoGradient
         } //end if phi
     } // end forall Points
 }

template<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)

todo{Whate are P1 and P2 over here?}

Todo:: this has recently been fixed

Definition at line 2294 of file StatisticsVector.hcc.

References PeriodicBoundary::getDistance(), PeriodicBoundary::shiftPositions(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

 {
     if (periodicWalls)
         for (unsigned int k = wp; k < boundaryHandler.getNumberOfObjects(); k++)
         {
             PeriodicBoundary* b0 = dynamic_cast<PeriodicBoundary*>(boundaryHandler.getObject(k));
             if (b0)
             {
                 b0->getDistance(P1);
                 b0->shiftPositions(P1, P2);
                 evaluate_force_statistics(k + 1);
                 b0->shiftPositions(P1, P2);
             }
         }
     
     //evaluate fields that only depend on CParticle parameters
     Mdouble phi; //Course graining integral
     for (unsigned int i = 0; i < Points.size(); i++)
     {
         phi = Points[i].CG_function(P);
         if (phi != 0.0)
         {
             Points[i].NormalTraction += P1_P2_NormalTraction * phi;
             Points[i].TangentialTraction += P1_P2_TangentialTraction * phi;
             if (getDoGradient())
             {
                 Vec3D d_phi = Points[i].CG_gradient(P, phi);
                 dx[i].NormalTraction += P1_P2_NormalTraction * d_phi.X;
                 dy[i].NormalTraction += P1_P2_NormalTraction * d_phi.Y;
                 dz[i].NormalTraction += P1_P2_NormalTraction * d_phi.Z;
                 dx[i].TangentialTraction += P1_P2_TangentialTraction * d_phi.X;
                 dy[i].TangentialTraction += P1_P2_TangentialTraction * d_phi.Y;
                 dz[i].TangentialTraction += P1_P2_TangentialTraction * d_phi.Z;
             } //end if getDoGradient
         } //end if phi
     } // end forall Points
     
 }

template<StatType T>

Mdouble StatisticsVector< T >::evaluateIntegral	(	Mdouble	n1,
		Mdouble	n2,
		Mdouble	t
	)

inline

Definition at line 908 of file StatisticsVector.h.

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

template<StatType T>

Mdouble StatisticsVector< T >::evaluatePolynomial ( Mdouble r )

inline

Definition at line 892 of file StatisticsVector.h.

     {
         return CGPolynomial.evaluate(r);
     }

template<StatType T>

Mdouble StatisticsVector< T >::evaluatePolynomialGradient ( Mdouble r )

inline

Definition at line 900 of file StatisticsVector.h.

     {
         return CGPolynomial.evaluateGradient(r);
     }

template<StatType T>

void StatisticsVector< T >::finishStatistics ( )

virtual

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

Reimplemented from DPMBase.

Definition at line 961 of file StatisticsVector.hcc.

 {
     if (getDoTimeAverage())
         write_time_average_statistics();
     statFile.close();
 }

template<StatType T>

void StatisticsVector< T >::gather_force_statistics_from_fstat_and_data ( )

get force statistics from particle collisions

Todo:: {fstat misses torque; any new implementation of storing forces should also store the torque}

Definition at line 2148 of file StatisticsVector.hcc.

References MULTIPLE_FILES, and MULTIPLE_FILES_PADDED.

 {
     //read in first three lines (should start with '#')
     std::string dummy;
     if (fStatFile.getFileType() == FileType::MULTIPLE_FILES || fStatFile.getFileType() == FileType::MULTIPLE_FILES_PADDED)
         fStatFile.openNextFile(std::fstream::in);
     getline(fStatFile.getFstream(), dummy);
     getline(fStatFile.getFstream(), dummy);
     getline(fStatFile.getFstream(), dummy);
     
     static Mdouble time;
     static int index1, index2;
     static Vec3D Contact;
     static Mdouble delta, ctheta, fdotn, fdott;
     static Vec3D P1_P2_normal_, P1_P2_tangential;
     static BaseParticle PI;
     static BaseParticle PJ;
     
     int counter = 0;
     //go through each line in the fstat file; break when eof or '#' or newline
     while ((fStatFile.getFstream().peek() != -1) && (fStatFile.getFstream().peek() != '#'))
     { //(!fStatFile.eof())
         counter++;
         /* # 1: time
          # 2: particle Number i
          # 3: contact partner j (particles >= 0, walls < 0)
          # 4: x-position \
          # 5: y-position  > of the contact point (I hope)
          # 6: z-position /
          # 7: delta = overlap at the contact
          # 8: ctheta = length of the tangential spring
          # 9: P1_P2_normal force |f^n|
          # 10: remaining (tangential) force |f^t|=|f-f^n|
          # 11-13: P1_P2_normal unit vector nx, ny, nz
          # 14-16: tangential unit vector tx, ty, tz
          */
         fStatFile.getFstream()
                 >> time
                 >> index1
                 >> index2
                 >> Contact
                 >> delta
                 >> ctheta
                 >> fdotn
                 >> fdott
                 >> P1_P2_normal_
                 >> P1_P2_tangential;
         fStatFile.getFstream().ignore(256, '\n');
         //Finished reading fstat file
         gatherContactStatistics(index1, index2, Contact, delta, ctheta, fdotn, fdott, P1_P2_normal_, P1_P2_tangential);
         
     }
     if (verbosity > 1)
         std::cout << "#forces=" << counter << ", t=" << time << std::endl;
 }

template<StatType T>

void StatisticsVector< T >::gather_force_statistics_from_p3c ( int version )

get force statistics from particle collisions

Todo:: tw: a better way to do p3/p4 statistics is to move p3/p4 files to Mercury files and perform statistics there

Todo:: This has to be tested

Definition at line 1522 of file StatisticsVector.hcc.

References Vec3D::dot(), Vec3D::getLength(), Vec3D::getLengthSquared(), BaseInteractable::getPosition(), BaseParticle::getRadius(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

 {
     std::string dummyStr;
     Mdouble N, dummy;
     
     //Ignoring first line
     getline(p3c_file, dummyStr);
     if (verbosity>2)
         std::cout << dummyStr << std::endl;
     //Reading second line
     p3c_file >> dummy;
     setTime(dummy);
     p3c_file >> N;
     std::cout << "t= " << getTime() << ": loading " << N << " contacts ..." << std::endl;
     getline(p3c_file, dummyStr);
     //Ignoring third line
     getline(p3c_file, dummyStr);
     if (verbosity>2)
         std::cout << dummyStr << std::endl;
     //Checking for end of file
     //if (p3c_file.eof()||p3c_file.peek()==-1) return false;
     
     //set up index list
     unsigned int Nmax = 0;
     for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
         Nmax = std::max(Nmax, (*it)->getIndex());
     static std::vector<int> index;
     index.resize(Nmax+1);
     {
         int i = 0;
         for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
         {
             index[(*it)->getIndex()] = i;
             i++;
         }
     }
     
     // ==> Particledata_4.2to4.25s.csv.p4c <==
     // TIMESTEP CONTACTS
     // 4.2001  31367
     //  P1 P2 CX CY CZ FX FY FZ
     //  9711 9424 0.0549396 0.009705380.149582 4.37283e-05 2.02347e-05 -0.000250662
     
     // ==> /Users/weinhartt/CarlosLabra/ConfinedCompression_P3P_Spheres.csv.p3c <==
     // TIMESTEP CONTACTS
     // 6.500010000000e-01  21886
     //  P1 P2 FX FY FZ
     //  1126 218 -1.899612000000e-02 6.716474700000e-02 -1.527340000000e-02
     
     int index1_p3, index2_p3;
     int index1, index2;
     Vec3D Contact, Force;
     Mdouble delta, fdotn, fdott;
     Vec3D P1_P2_normal, P1_P2_tangential;
     Mdouble dist;
     BaseParticle* P1;
     BaseParticle* P2;
     
     Vec3D contact;
     int counter = 0;
     //go through N lines in the p4c file
     for (int i = 0; i < N; i++)
     {
         counter++;
         if (version == 3)
         {
             p3c_file
                     >> index1_p3
                     >> index2_p3
                     >> Force;
         }
         else //p4p
         {
             p3c_file
                     >> index1_p3
                     >> index2_p3
                     >> Contact
                     >> Force;
         }
         index1 = index[index1_p3];
         index2 = index[index2_p3];
         P1 = particleHandler.getObject(index1);
         P2 = particleHandler.getObject(index2);
         P1_P2_normal = P1->getPosition() - P2->getPosition();
         // deal with periodic boundaries; correct if the distance between two particles is bigger than half the domain size
         // (only works if periodic walls are >2d apart)
         std::cout << "#forces="  << std::setprecision(9) << P1_P2_normal.Y << "|" << getYMax()-getYMin() << std::endl;
         if (2.0*P1_P2_normal.X > getXMax()-getXMin())
             P1_P2_normal.X -=     getXMax()-getXMin();
         else if (2.0*P1_P2_normal.X < getXMin()-getXMax())
             P1_P2_normal.X +=     getXMax()-getXMin();
 
         else if (2.0*P1_P2_normal.Y > getYMax()-getYMin())
             P1_P2_normal.Y -=     getYMax()-getYMin();
         else if (2.0*P1_P2_normal.Y < getYMin()-getYMax())
             P1_P2_normal.Y +=     getYMax()-getYMin();
 
         else if (2.0*P1_P2_normal.Z > getZMax()-getZMin())
             P1_P2_normal.Z -=     getZMax()-getZMin();
         else if (2.0*P1_P2_normal.Z < getZMin()-getZMax())
             P1_P2_normal.Z +=     getZMax()-getZMin();
         std::cout << "#forces="  << std::setprecision(9) << P1_P2_normal.Y << "|" << getYMax()-getYMin() << std::endl;
 
         dist = Vec3D::getLength(P1_P2_normal);
         P1_P2_normal /= dist;
         delta = P1->getRadius() + P2->getRadius() - dist;
         if (version == 3)
         {
             Contact = P2->getPosition() + P1_P2_normal * (P2->getRadius() - delta / 2);
         }
         fdotn = -Vec3D::dot(Force, P1_P2_normal) / Vec3D::getLengthSquared(P1_P2_normal);
         P1_P2_tangential = Force + fdotn * P1_P2_normal; //get tangential force
         fdott = Vec3D::getLength(P1_P2_tangential);
         if (fdott == 0)
         { 
             P1_P2_tangential = Vec3D(0, 0, 0);
         }
         else
         {
             P1_P2_tangential /= fdott;
         }
         //Finished reading file
         std::cout << "#forces=" << P1_P2_normal << "|" << fdotn << std::endl;
         gatherContactStatistics(index1, index2, Contact, delta, 0, fdotn, fdott, P1_P2_normal, -P1_P2_tangential);
         gatherContactStatistics(index2, index1, Contact, delta, 0, fdotn, fdott, -P1_P2_normal, P1_P2_tangential);
     }
     
     if (verbosity > 1)
         std::cout << "#forces=" << counter << std::endl;
     //Ignoring last line
     getline(p3c_file, dummyStr);
     if (verbosity>2)
         std::cout << dummyStr << std::endl;
     
     gather_force_statistics_from_p3w(version, index);
 }

template<StatType T>

void StatisticsVector< T >::gather_force_statistics_from_p3w	(	int	version,
		std::vector< int > &	index
	)

get force statistics from particle collisions

Todo:: This has to be tested

Definition at line 1662 of file StatisticsVector.hcc.

References Vec3D::dot(), Vec3D::getLength(), Vec3D::getLengthSquared(), BaseInteractable::getPosition(), and BaseParticle::getRadius().

 {
     std::string dummyStr;
     Mdouble N, dummy;
     
     //Ignoring first line
     getline(p3w_file, dummyStr);
     if (verbosity>2)
         std::cout << dummyStr << std::endl;
     //Reading second line
     p3w_file >> dummy;
     setTime(dummy);
     p3w_file >> N;
     std::cout << "t= " << getTime() << ": loading " << N << " wall contacts ..." << std::endl;
     getline(p3w_file, dummyStr);
     //Ignoring third line
     getline(p3w_file, dummyStr);
     if (verbosity>2)
         std::cout << dummyStr << std::endl;
     
     // ==> /Users/weinhartt/CarlosLabra/ConfinedCompression_P3P_Spheres.csv.p3w <==
     // TIMESTEP CONTACTS
     // 6.500010000000e-01  3778
     //  P1 FX FY FZ NX NY NZ (N=Branch)
     //  671 0.000000000000e+00 0.000000000000e+00 1.965970000000e-01 0.000000000000e+00 0.000000000000e+00 -1.991510000000e-03
     //  ==> Particledata_4.2to4.25s.csv.p4w <==
     // TIMESTEP CONTACTS
     // 4.2001  5770
     //  P1 CX CY CZ FX FY FZ
     //   10590 0.00944221 0.0120.130307 -3.14611e-06 -0.000512936 0.000102539
     //   10146 0.0615829 00.162907 2.02968e-05 0.000948522 0.000188616
     
     int index1_p3;
     int index1, index2;
     Vec3D Contact, Force, Branch;
     Mdouble delta, fdotn, fdott;
     Vec3D P1_P2_normal, P1_P2_tangential;
     Mdouble dist;
     BaseParticle* P1;
     
     Vec3D contact;
     int counter = 0;
     //go through each line in the fstat file; break when eof or '#' or newline
     for (int i = 0; i < N; i++)
     {
         counter++;
         if (version == 3)
         {
             p3w_file
                     >> index1_p3
                     >> Force
                     >> Branch;
         }
         else
         {
             p3w_file
                     >> index1_p3
                     >> Contact
                     >> Force;
         }
         index1 = index[index1_p3];
         index2 = -1;
         P1 = particleHandler.getObject(index1);
         if (version == 3)
         {
             //Branch *= -1;
             Contact = P1->getPosition() + Branch;
         }
         else
         {
             Branch = Contact - P1->getPosition();
         }
         dist = Vec3D::getLength(Branch);
         P1_P2_normal = Branch / dist;
         delta = P1->getRadius() - dist; //check
         fdotn = -Vec3D::dot(Force, P1_P2_normal) / Vec3D::getLengthSquared(P1_P2_normal);
         P1_P2_tangential = Force + fdotn * P1_P2_normal; //get tangential force
         fdott = Vec3D::getLength(P1_P2_tangential);
         if (fdott == 0)
         { 
             P1_P2_tangential = Vec3D(0, 0, 0);
         }
         else
         {
             P1_P2_tangential /= fdott;
         }
         //Finished reading fstat file
         gatherContactStatistics(index1, index2, Contact, delta, 0, fdotn, fdott, P1_P2_normal, -P1_P2_tangential);
     }
     
     if (verbosity > 1)
         std::cout << "#wall forces=" << counter << std::endl;
     //Ignoring last line
     getline(p3w_file, dummyStr);
     if (verbosity>2)
         std::cout << dummyStr << std::endl;
 }

template<StatType T>

void StatisticsVector< T >::gatherContactStatistics	(	unsigned int index1	,
		int index2	,
		Vec3D Contact	,
		Mdouble delta	,
		Mdouble ctheta	UNUSED,
		Mdouble fdotn	,
		Mdouble fdott	,
		Vec3D P1_P2_normal_	,
		Vec3D P1_P2_tangential
	)

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 such that it always extends to z->inf

2012Nov28TW I changed this from P1 to P2; was this mistake always there?

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 DPMBase.

Definition at line 1985 of file StatisticsVector.hcc.

References Vec3D::dot(), Matrix3D::dyadic(), MatrixSymmetric3D::selfDyadic(), and Vec3D::Z.

 {
     Vec3D P1_P2_VelocityAverage, P1_P2_VelocityDifference, P1_P2_Force;
     Mdouble norm_dist;
     if (check_current_time_for_statistics())
     {
         //make sure that for a particle-fixed particle collision the fixed particle is on index2 (switch particles if needed)
         if (particleHandler.getObject(index1)->isFixed())
         {
             int tmp = index1;
             index1 = index2;
             index2 = tmp;
             P1_P2_normal_ *= -1;
         }
         
         if (satisfiesInclusionCriteria(particleHandler.getObject(index1)))
         {
             P1_P2_normal = P1_P2_normal_;
             
             if (getSystemDimensions() == 2)
                 Contact.Z = 0.0;
             
             if (index2 < 0)
             { //wall-particle collision
                 P1_P2_distance = particleHandler.getObject(index1)->getRadius() - delta;
                 P2 = Contact;
                 P1 = Contact - P1_P2_normal * P1_P2_distance; //note, here we use a minus
                 norm_dist = P1_P2_distance;
                 P1_P2_VelocityAverage = (particleHandler.getObject(index1)->getVelocity()) / 2;
                 P1_P2_VelocityDifference = (particleHandler.getObject(index1)->getVelocity());
                 
                 if (StressTypeForFixedParticles == 3)
                 {
                     P1_P2_distance += 300;
                     P2 += 300 * P1_P2_normal;
                     //this is because wall-particle collisions appear only once in fstat
                     norm_dist = P1_P2_distance;
                 }
             }
             else if (particleHandler.getObject(index2)->isFixed() || !satisfiesInclusionCriteria(particleHandler.getObject(index2)))
             { //particle-fixed particle collision (external force acts at contact point)
                 if (particleHandler.getObject(index2)->isFixed() && StressTypeForFixedParticles == 3)
                 {
                     //infinite extension of stress
                     P1_P2_distance = particleHandler.getObject(index1)->getRadius() + particleHandler.getObject(index2)->getRadius() - delta;
                     P1 = Contact + P1_P2_normal * (0.5 * P1_P2_distance); //flowing particle
                     P2 = Contact - P1_P2_normal * (0.5 * P1_P2_distance); //fixed particle
                     if (P1_P2_normal.Z < 0)
                     { //in this case revert flowing and fixed particle
                       //std::cerr << "Warning: Force normal upwards on base z=" << P1.Z << "." << P2.Z << std::endl;
                       //turning force the other way to achieve sigma=0 at z=inf
                         P1_P2_normal *= -1.0;
                         fdotn *= -1.0;
                     }
                     P1_P2_distance = setInfinitelyLongDistance();
                     P2 = P1 - P1_P2_normal * P1_P2_distance; //move endpoint of force line downwards beyond fixed particles out of the domain
                     if (P1_P2_normal.Z < 0)
                         std::cout << P1_P2_distance << " z" << P1.Z << " " << P2.Z << " f" << fdotn << std::endl;
                 }
                 else if (StressTypeForFixedParticles == 1 || (!particleHandler.getObject(index2)->isFixed() && StressTypeForFixedParticles == 3))
                 {
                     //add force from flowing particle to contact point to stress
                     P1_P2_distance = particleHandler.getObject(index1)->getRadius() + particleHandler.getObject(index2)->getRadius() - delta;
                     P1 = Contact + P1_P2_normal * (0.5 * P1_P2_distance); //1st particle
                     //P2 = Contact; //Contact point;
                     //P1_P2_distance /= 2.0;
                     P1_P2_distance = particleHandler.getObject(index1)->getRadius() - delta / 2; //Contact point;
                     P2 = P1 - P1_P2_normal * P1_P2_distance; //Contact point (corrected for polydispersed particles);
                 }
                 else
                 {
                     //add force from flowing particle to fixed particle to stress
                     P1_P2_distance = particleHandler.getObject(index1)->getRadius() + particleHandler.getObject(index2)->getRadius() - delta;
                     P1 = Contact + P1_P2_normal * (0.5 * P1_P2_distance);
                     P2 = Contact - P1_P2_normal * (0.5 * P1_P2_distance);
                 }
                 norm_dist = P1_P2_distance; //*2.0*Particles[index1].Radius/(Particles[index2].Radius+Particles[index1].Radius);
                 //This is the length of the contact in particle 1.
                 P1_P2_VelocityAverage = particleHandler.getObject(index2)->getVelocity() / 2;
                 P1_P2_VelocityDifference = -particleHandler.getObject(index2)->getVelocity();
             }
             else
             { //particle-particle collision
             
                 if (particleHandler.getObject(index2)->isFixed())
                     std::cout << "ERROR" << std::endl;
                 P1_P2_distance = particleHandler.getObject(index1)->getRadius() + particleHandler.getObject(index2)->getRadius() - delta;
                 
                 P1 = Contact + P1_P2_normal * (0.5 * P1_P2_distance);
                 P2 = Contact - P1_P2_normal * (0.5 * P1_P2_distance);
                 //This is the length of the contact in particle 1.
                 norm_dist = P1_P2_distance * particleHandler.getObject(index1)->getRadius() / (particleHandler.getObject(index2)->getRadius() + particleHandler.getObject(index1)->getRadius());
                 P1_P2_VelocityAverage = (particleHandler.getObject(index1)->getVelocity() + particleHandler.getObject(index2)->getVelocity()) / 2;
                 P1_P2_VelocityDifference = (particleHandler.getObject(index1)->getVelocity() - particleHandler.getObject(index2)->getVelocity());
             }
             //Particles[max(index1,index2)] is chosen so that it is the non-wall, non-fixed particle
             P1_P2_Fabric = MatrixSymmetric3D::selfDyadic(P1_P2_normal) * particleHandler.getObject(index1)->getVolume();
             //fix
             //~ P1_P2_Fabric = SymmetrizedDyadic(P1_P2_normal, P1_P2_normal) * particleHandler.getObject(max(index1,index2))->getVolume(speciesHandler.getObject());
             P1_P2_NormalStress = Matrix3D::dyadic(P1_P2_normal, P1_P2_normal) * (fdotn * norm_dist);
             P1_P2_TangentialStress = Matrix3D::dyadic(P1_P2_tangential, P1_P2_normal) * (fdott * norm_dist);
             P1_P2_NormalTraction = P1_P2_normal * fdotn;
             P1_P2_TangentialTraction = P1_P2_tangential * fdott;
             P1_P2_ContactCoupleStress = Matrix3D::dyadic(P1_P2_NormalTraction + P1_P2_TangentialTraction, -0.5 * P1_P2_normal * norm_dist);
             P1_P2_Contact = Contact;
             //P1_P2_Potential = (speciesHandler.getObject(0)->getStiffness() * mathsFunc::square(delta) + dynamic_cast<FrictionalSpecies*>(speciesHandler.getObject(0))->getSlidingStiffness() * mathsFunc::square(ctheta)) / 2;
             P1_P2_Potential = 0.0;
             P1_P2_Force = fdotn * P1_P2_normal + fdott * P1_P2_tangential;
             P1_P2_Dissipation = Vec3D::dot(P1_P2_VelocityDifference, P1_P2_Force) / 2;
             P1_P2_CollisionalHeatFlux = P1_P2_normal * (P1_P2_distance / 2 * Vec3D::dot(P1_P2_VelocityAverage, fdotn * P1_P2_normal + fdott * P1_P2_tangential)) + P1_P2_Potential * P1_P2_VelocityAverage;
             //Now P1_P2 distance and normal are in non-averaged directions - This is the projection of the distance on to the CG directions. Required for the later statistics hence done now, not before.
             P1_P2_distance *= sqrt(Points[0].dotNonAveraged(P1_P2_normal, P1_P2_normal));
             
             //If the normal to the contact and the averaging direction are EXACTLY+- parallel then set the normal to zero. Otherwise you would end up with nan, which is not good.
             if (P1_P2_distance > 1e-12)
                 P1_P2_normal = (P1 - P2) / P1_P2_distance;
             else
                 P1_P2_normal = Vec3D(0.0, 0.0, 0.0);
             //P1_P2_normal = (P1-P2)/P1_P2_distance;
             
             //if 2nd particle is wall, fixed particle or not-included particle
             if (index2 < 0 || particleHandler.getObject(index2)->isFixed() || !satisfiesInclusionCriteria(particleHandler.getObject(index2)))
             {
                 // << "/" << P2 << std::endl;
                 if (StressTypeForFixedParticles != 0)
                 { //only if wall - flow contact adds anything to the stress
                     evaluate_force_statistics();
                 }
                 if (StressTypeForFixedParticles != 3 || !(index2 < 0 || particleHandler.getObject(index2)->isFixed()))
                 { //only if stress is not infinitely extended
                     Vec3D P;
                     if (StressTypeForFixedParticles == 0)
                     {
                         P = P1; //Traction at flow particle
                     }
                     else
                     {
                         P = P2; //Traction at contact point, resp wall particle
                     }
                     evaluate_wall_force_statistics(P);
                 }
             }
             else
             {
                 evaluate_force_statistics();
             }
         }
     }
 }

template<StatType T>

std::vector<StatisticsPoint<T> > StatisticsVector< T >::getCGPoints ( )

inline

Definition at line 672 of file StatisticsVector.h.

     {
         return Points;
     }

template<StatType T>

CG StatisticsVector< T >::getCGShape ( )

inline

Definition at line 288 of file StatisticsVector.h.

     {
         return CG_type;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getCGTimeAveragingInterval ( )

inline

Definition at line 261 of file StatisticsVector.h.

     {
         return tintStat;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getCGTimeMin ( )

inline

Definition at line 242 of file StatisticsVector.h.

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

     {
         return tminStat;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getCGWidth ( )

inline

Definition at line 326 of file StatisticsVector.h.

     {
         return std::sqrt(w2);
     }

template<StatType T>

Mdouble StatisticsVector< T >::getCGWidth_over_rmax ( )

inline

Definition at line 583 of file StatisticsVector.h.

     {
         return w_over_rmax;
     }

template<StatType T>

bool StatisticsVector< T >::getCGWidthalls ( )

inline

Definition at line 550 of file StatisticsVector.h.

     {
         return walls;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getCGWidthSquared ( )

inline

Definition at line 334 of file StatisticsVector.h.

     {
         return w2;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getCutoff ( )

inline

Definition at line 342 of file StatisticsVector.h.

     {
         return cutoff;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getCutoff2 ( )

inline

Definition at line 350 of file StatisticsVector.h.

     {
         return mathsFunc::square(cutoff);
     }

template<StatType T>

bool StatisticsVector< T >::getDoDoublePoints ( )

inline

Definition at line 844 of file StatisticsVector.h.

     {
         return doDoublePoints;
     }

template<StatType T>

bool StatisticsVector< T >::getDoGradient ( )

inline

Definition at line 478 of file StatisticsVector.h.

     {
         return doGradient;
     }

template<StatType T>

bool StatisticsVector< T >::getDoIgnoreFixedParticles ( )

inline

Definition at line 510 of file StatisticsVector.h.

     {
         return ignoreFixedParticles;
     }

template<StatType T>

bool StatisticsVector< T >::getDoPeriodicWalls ( )

inline

Definition at line 566 of file StatisticsVector.h.

     {
         return periodicWalls;
     }

template<StatType T>

bool StatisticsVector< T >::getDoTimeAverage ( )

inline

Definition at line 400 of file StatisticsVector.h.

     {
         return doTimeAverage;
     }

template<StatType T>

bool StatisticsVector< T >::getDoVariance ( )

inline

Definition at line 461 of file StatisticsVector.h.

     {
         return doVariance;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getMirrorAtDomainBoundary ( )

inline

Definition at line 444 of file StatisticsVector.h.

     {
         return mirrorAtDomainBoundary;
     }

template<StatType T>

void StatisticsVector< T >::getN	(	int &	nx_,
		int &	ny_,
		int &	nz_
	)

inline

Definition at line 303 of file StatisticsVector.h.

     {
         nx_ = nx;
         ny_ = ny;
         nz_ = nz;
     }

template<StatType T>

int StatisticsVector< T >::getNTimeAverage ( )

inline

Definition at line 794 of file StatisticsVector.h.

     {
         return nTimeAverage;
     }

template<StatType T>

bool StatisticsVector< T >::getNTimeAverageReset ( )

inline

Definition at line 860 of file StatisticsVector.h.

     {
         return nTimeAverageReset;
     }

template<StatType T>

int StatisticsVector< T >::getNX ( )

inline

Definition at line 194 of file StatisticsVector.h.

     {
         return nx;
     }

template<StatType T>

int StatisticsVector< T >::getNY ( )

inline

Definition at line 202 of file StatisticsVector.h.

     {
         return ny;
     }

template<StatType T>

int StatisticsVector< T >::getNZ ( )

inline

Definition at line 210 of file StatisticsVector.h.

     {
         return nz;
     }

template<StatType T>

std::string StatisticsVector< T >::getPolynomialName ( )

inline

Definition at line 828 of file StatisticsVector.h.

     {
         return CGPolynomial.getName();
     }

template<StatType T>

unsigned int StatisticsVector< T >::getStepSize ( ) const

inline

Definition at line 916 of file StatisticsVector.h.

     {
         return stepSize_;
     }

template<StatType T>

int StatisticsVector< T >::getStressTypeForFixedParticles ( )

inline

Definition at line 417 of file StatisticsVector.h.

     {
         return StressTypeForFixedParticles;
     }

template<StatType T>

bool StatisticsVector< T >::getSuperExact ( )

inline

Definition at line 494 of file StatisticsVector.h.

     {
         return superexact;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getTimeMaxStat ( )

inline

Definition at line 250 of file StatisticsVector.h.

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

     {
         if (std::isnan(tmaxStat))
             return getTimeMax();
         else
             return tmaxStat;
     }

template<StatType T>

Vec3D StatisticsVector< T >::getVelocityProfile ( Vec3D Position )

protected

Definition at line 366 of file StatisticsVector.hcc.

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

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

template<StatType T>

int StatisticsVector< T >::getVerbosityLevel ( ) const

inline

Definition at line 534 of file StatisticsVector.h.

     {
         return verbosity;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getXMaxStat ( )

inline

Definition at line 713 of file StatisticsVector.h.

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

     {
         if (std::isnan(xmaxStat))
             return getXMax();
         else
             return xmaxStat;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getXMinStat ( )

inline

Functions to acces and change the domain of statistics.

Definition at line 680 of file StatisticsVector.h.

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

     {
         if (std::isnan(xminStat))
             return getXMin();
         else
             return xminStat;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getYMaxStat ( )

inline

Definition at line 724 of file StatisticsVector.h.

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

     {
         if (std::isnan(ymaxStat))
             return getYMax();
         else
             return ymaxStat;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getYMinStat ( )

inline

Definition at line 691 of file StatisticsVector.h.

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

     {
         if (std::isnan(yminStat))
             return getYMin();
         else
             return yminStat;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getZMaxStat ( )

inline

Definition at line 735 of file StatisticsVector.h.

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

     {
         if (std::isnan(zmaxStat))
             return getZMax();
         else
             return zmaxStat;
     }

template<StatType T>

Mdouble StatisticsVector< T >::getZMinStat ( )

inline

Definition at line 702 of file StatisticsVector.h.

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

     {
         if (std::isnan(zminStat))
             return getZMin();
         else
             return zminStat;
     }

template<StatType T>

void StatisticsVector< T >::initialiseStatistics ( )

virtual

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

Reimplemented from DPMBase.

Definition at line 883 of file StatisticsVector.hcc.

References StatisticsVector< T >::setCGWidth2(), and StatisticsVector< T >::setPositions().

 {
     reset_statistics();
     
     StatisticsVector<T>::setCGWidth2(StatisticsVector<T>::getCGWidthSquared());
     
     StatisticsVector<T>::setPositions();
     
     //set tmin and tmax if tint is set
     if (getCGTimeAveragingInterval() != 0.0)
     {
         setCGTimeMin(getTime());
         setTimeMaxStat(getCGTimeMin() + getCGTimeAveragingInterval());
         std::cout << "tint set, time goes from: " << getCGTimeMin() << " to: " << getTimeMaxStat() << std::endl;
     }
 
     statFile.open();
 
     statFile.getFstream() << Points.begin()->write_variable_names() << std::endl;
     statFile.getFstream() << print_CG() << std::endl;
     
     //std::couts variables
     if (verbosity > 1)
         std::cout << std::endl;
     if (verbosity > 0)
         std::cout << std::endl << printStat() << std::endl;
     
 }

template<StatType T>

void StatisticsVector< T >::jump_fstat ( )

Todo:: make sure the counters for fstat and data are used right here (check for consistency)

Definition at line 1953 of file StatisticsVector.hcc.

References MULTIPLE_FILES, and MULTIPLE_FILES_PADDED.

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

template<StatType T>

void StatisticsVector< T >::jump_p3c ( )

get force statistics from particle collisions

Definition at line 1481 of file StatisticsVector.hcc.

 {
     std::string dummyStr;
     Mdouble N, dummy;
     
     //Ignoring first line
     getline(p3c_file, dummyStr);
     std::cout << dummyStr << std::endl;
     //Reading second line
     p3c_file >> dummy;
     setTime(dummy);
     p3c_file >> N;
     std::cout << "ignoring t= " << getTime() << ": loading " << N << " contacts ..." << std::endl;
     getline(p3c_file, dummyStr);
     //Ignoring third line
     getline(p3c_file, dummyStr);
     std::cout << dummyStr << std::endl;
     //go through N lines in the p3c file
     for (int i = 0; i < N; i++)
         getline(p3c_file, dummyStr);
     
     //Ignoring first line
     getline(p3w_file, dummyStr);
     if (verbosity>2)
         std::cout << dummyStr << std::endl;
     //Reading second line
     p3w_file >> dummy;
     setTime(dummy);
     p3w_file >> N;
     std::cout << "ignoring t= " << getTime() << ": loading " << N << " wall contacts ..." << std::endl;
     getline(p3w_file, dummyStr);
     //Ignoring third line
     getline(p3w_file, dummyStr);
     //go through N lines in the p3w file
     for (int i = 0; i < N; i++)
         getline(p3w_file, dummyStr);
     
 }

template<StatType T>

bool StatisticsVector< T >::loadPositions ( const char * filename )

protected

Definition at line 1927 of file StatisticsVector.hcc.

 {
     //This opens the file the data will be recalled from
     std::fstream is;
     is.open(filename, std::fstream::in);
     if (!is.is_open() || is.bad())
     {
         std::cout << "Loading positions file " << filename << " failed" << std::endl;
         return false;
     }
 
     double N=0;
     is >> N;
     positions_.resize(N);
     for (int i=0; i<N; i++)
     {
         is >> positions_[i];
     }
     //Close the file
     is.close();
     std::cout << "Loaded " << positions_.size() << " positions from file " << filename << std::endl;
     return true;
 }

template<StatType T>

bool StatisticsVector< T >::loadVelocityProfile ( const char * filename )

protected

Definition at line 304 of file StatisticsVector.hcc.

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

 {
     //This opens the file the data will be recalled from
     std::fstream is;
     is.open(filename, std::fstream::in);
     if (!is.is_open() || is.bad())
     {
         std::cout << "Loading data file " << filename << " failed" << std::endl;
         return false;
     }
     
     //Read the file; we assume that the stat file we read in has the same stattype and domain size as the problem.
     std::string line_string;
     getline(is, line_string);
     std::string dummy;
     is >> dummy >> dummy;
     is >> dummy >> dummy;
     double xmin, xmax, ymin, ymax, zmin, zmax;
     is >> dummy >> xmin >> xmax >> ymin >> ymax >> zmin >> zmax;
     VelocityProfile_Min = Vec3D(xmin, ymin, zmin);
     int nx, ny, nz;
     is >> dummy >> nx >> ny >> nz;
     double dx = (xmax - xmin) / (nx - 1);
     double dy = (ymax - ymin) / (ny - 1);
     double dz = (zmax - zmin) / (nz - 1);
     VelocityProfile_D = Vec3D(dx, dy, dz);
     std::string statType;
     is >> dummy >> statType;
     getline(is, line_string);
     getline(is, line_string);
     
     setNX(nx);
     setNY(ny);
     setNZ(nz);
     setXMinStat(xmin);
     setYMinStat(ymin);
     setZMinStat(zmin);
     setXMaxStat(xmax);
     setYMaxStat(ymax);
     setZMaxStat(zmax);
     
     VelocityProfile.resize(nx * ny * nz);
     Vec3D M;
     double rho;
     int n = 0;
     for (int i = 0; i < nx; i++)
         for (int j = 0; j < ny; j++)
             for (int k = 0; k < nz; k++)
             {
                 is >> dummy >> dummy >> dummy >> dummy >> rho >> M.X >> M.Y >> M.Z;
                 getline(is, line_string);
                 VelocityProfile[n] = M / rho;
                 n++;
             }
     
     //Close the file
     is.close();
     std::cout << "Loaded data file " << filename << std::endl;
     return true;
 }

template<StatType T>

void StatisticsVector< T >::outputStatistics ( )

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 DPMBase.

Definition at line 2336 of file StatisticsVector.hcc.

 {
     if (check_current_time_for_statistics())
     {
         for (unsigned int i = 0; i < Points.size(); i++)
             Points[i].setCGInverseVolume();
         
         for (std::vector<BaseParticle*>::iterator P = particleHandler.begin(); P != particleHandler.end(); ++P)
         {
             //unfix particles (turn off to ignore fixed particles)
             if ((!ignoreFixedParticles) && (*P)->isFixed())
             {
                 (*P)->unfix();
             }
             //ignore fixed particles
             
             if (satisfiesInclusionCriteria(*P) && !(*P)->isFixed())
             {
                 evaluate_particle_statistics(P);
             }
             
         }
     }
 }

template<StatType T>

std::string StatisticsVector< T >::print_CG ( )

Output coarse graining variables.

Output names of statistical variables.

Todo:: {Thomas: stat file should include all variables such as CGShape Stattype gradient timevariance timeaveraged, ...}

Definition at line 827 of file StatisticsVector.hcc.

References Polynomial.

 {
     std::stringstream ss;
     ss << "w " << sqrt(getCGWidthSquared())
             << " dim " << getSystemDimensions()
             << " domainStat " << getXMinStat() << " " << getXMaxStat()
             << " " << getYMinStat() << " " << getYMaxStat()
             << " " << getZMinStat() << " " << getZMaxStat()
             << " n " << nx << " " << ny << " " << nz
             << " statType " << statType;
     if (CG_type == Polynomial)
         ss << " CGPolynomial " << CGPolynomial;
     else
         ss << " CG_type " << CG_type << " cutoff " << cutoff;
     if (doVariance)
         ss << " doVariance";
     if (doGradient)
         ss << " doGradient";
     if (doTimeAverage)
     {
         ss << " doTimeAverage";
         if (nTimeAverageReset != -1)
             ss << " nTimeAverageReset " << nTimeAverageReset;
     }
     if (indSpecies != -1)
         ss << " indSpecies " << indSpecies;
     if (rmin != 0)
         ss << " rmin " << rmin;
     if (rmax != 1e20)
         ss << " rmax " << rmax;
     if (hmax != 1e20)
         ss << " hmax " << hmax;
     if (!walls)
         ss << " noWalls";
     if (!periodicWalls)
         ss << " noPeriodicWalls";
     if (!ignoreFixedParticles)
         ss << " includeFixedParticles";
     if (StressTypeForFixedParticles != 1)
         ss << " StressTypeForFixedParticles " << StressTypeForFixedParticles;
     return ss.str();
 }

template<StatType T>

void StatisticsVector< T >::print_help ( )

Definition at line 654 of file StatisticsVector.hcc.

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

template<StatType T>

std::string StatisticsVector< T >::printStat ( )

Outputs member variable values to a std::string.

Outputs StatisticsVector.

Definition at line 781 of file StatisticsVector.hcc.

References Gaussian, HeavisideSphere, and Polynomial.

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

template<StatType T>

void StatisticsVector< T >::processStatistics ( bool usethese )

virtual

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

Reimplemented from DPMBase.

Definition at line 969 of file StatisticsVector.hcc.

 {
     if (check_current_time_for_statistics() && usethese)
     {
         if (getMirrorAtDomainBoundary() != 0.0)
         {
             for (unsigned int i = 0; i < Points.size(); i++)
             {
                 if (Points[i].mirrorParticle >= 0)
                 {
                     //add values to the mirrorParticle
                     Points[Points[i].mirrorParticle] += Points[i];
                     Points[i].set_zero();
                     if (getDoGradient())
                     {
                         dx[Points[i].mirrorParticle] += dx[i];
                         dy[Points[i].mirrorParticle] += dy[i];
                         dz[Points[i].mirrorParticle] += dz[i];
                     }
                 }
             }
         }
         //~ for (unsigned int i=0; i<Points.size(); i++) 
         //~ Points[i].Displacement /= mathsFunc::square(Points[i].Density);
         if (getDoTimeAverage())
         {
             for (unsigned int i = 0; i < timeAverage.size(); i++)
             {
                 timeAverage[i] += Points[i];
                 if (getDoVariance())
                     timeVariance[i] += Points[i].getSquared();
                 if (getDoGradient())
                 {
                     dxTimeAverage[i] += dx[i];
                     dyTimeAverage[i] += dy[i];
                     dzTimeAverage[i] += dz[i];
                 }
             }
             nTimeAverage++;
             if (nTimeAverage == nTimeAverageReset)
             {
                 std::cout << "write time average" << std::endl;
                 write_time_average_statistics();
                 nTimeAverage = 0;
                 for (unsigned int i = 0; i < timeAverage.size(); i++)
                     timeAverage[i].set_zero();
             }
         }
         else
         {
             write_statistics();
         }
     }
     reset_statistics();
 }

template<StatType T>

bool StatisticsVector< T >::read_next_from_p3p_file ( )

protected

Definition at line 1266 of file StatisticsVector.hcc.

References constants::pi.

 {
     //uncomment if you want data file output
     //MD::outputXBallsData();
     
     unsigned int N;
     Mdouble dummy;
     std::string dummyStr;
     
     //read first parameters and check if we reached the end of the file
     
     //Ignoring first line
     getline(p3p_file, dummyStr);
     if (verbosity>2)
         std::cout << dummyStr << std::endl;
     //Reading second line
     p3p_file >> dummy;
     setTime(dummy);
     p3p_file >> N;
     std::cout << "t= " << getTime() << ": loading " << N << " particles ..." << std::endl;
     getline(p3p_file, dummyStr);
     //std::cout << dummyStr << std::endl;
     //Ignoring third line
     getline(p3p_file, dummyStr);
     if (verbosity>2)
         std::cout << dummyStr << std::endl;
     if (p3p_file.eof() || p3p_file.peek() == -1)
     {
         std::cout << "reached end of p3p file" << std::endl;
         return false;
     }
     
     BaseParticle P0;
     if (particleHandler.getNumberOfObjects() < N)
         while (particleHandler.getNumberOfObjects() < N)
             particleHandler.copyAndAddObject(P0);
     else
         while (particleHandler.getNumberOfObjects() > N)
             particleHandler.removeLastObject();
     
     //Read forth to (N+3)rd line
     Vec3D dummyVec;
     Vec3D position, velocity, angle, angularVelocity;
     int dummyInt;
     for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
     {
         p3p_file >> dummyInt;
         (*it)->setIndex(dummyInt); //ID
         //std::cout << "  ID " << (*it)->getIndex() << std::endl;
         p3p_file >> dummyInt;
         (*it)->setIndSpecies(dummyInt - 1); //Group
         p3p_file >> dummy;
         (*it)->setRadius(pow(dummy / (4. / 3. * constants::pi), 1. / 3.)); //Volume
         p3p_file >> dummy;
         (*it)->setMass(dummy); //Mass
         p3p_file
                 >> position
                 >> velocity
                 >> angularVelocity;
         (*it)->setPosition(position);
         (*it)->setVelocity(velocity);
         (*it)->setAngularVelocity(angularVelocity);
         //clean up tangential springs
     } //end for all particles
       //auto_setSystemDimensions();
     
     getline(p3p_file, dummyStr);
     if (verbosity>2)
         std::cout << dummyStr << std::endl;
     
     //fix particles, if data_FixedParticles!=0
     
     return true;
 }

template<StatType T>

bool StatisticsVector< T >::readNextDataFile ( unsigned int format )

Definition at line 913 of file StatisticsVector.hcc.

References DPMBase::readNextDataFile().

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

template<StatType T>

void StatisticsVector< T >::readStatArguments	(	int	argc,
		char *	argv[]
	)

Todo:: this should also set the ene,restart filename

Definition at line 398 of file StatisticsVector.hcc.

References Gaussian, and HeavisideSphere.

 {
     //check for options (standardly '-option argument')
     for (int i = 2; i < argc; i += 2)
     {
         std::cout << "interpreting input argument " << argv[i];
         if (i + 1 < argc)
             std::cout << " " << argv[i + 1];
         std::cout << std::endl;
         
         if (!strcmp(argv[i], "-CGtype"))
         {
             //CG_type_todo
             if (!strcmp(argv[i + 1], "HeavisideSphere"))
             {
                 setCGShape(HeavisideSphere);
             }
             else if (!strcmp(argv[i + 1], "Gaussian"))
             {
                 setCGShape(Gaussian);
             }
             else
             {
                 setCGShape(argv[i + 1]);
             }
         }
         else if (!strcmp(argv[i], "-w"))
         {
             double old = getCGWidth();
             setCGWidth(atof(argv[i + 1]));
             std::cout << " w changed from " << old << " to " << getCGWidth() << std::endl;
         }
         else if (!strcmp(argv[i], "-help"))
         {
             print_help();
             i--; //requires no argument
         }
         else if (!strcmp(argv[i], "-velocityprofile"))
         {
             loadVelocityProfile(argv[i + 1]);
         }
         else if (!strcmp(argv[i], "-h"))
         {
             set_h(atof(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-hx"))
         {
             set_hx(atof(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-hy"))
         {
             set_hy(atof(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-hz"))
         {
             set_hz(atof(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-n"))
         {
             setN(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-nx"))
         {
             setNX(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-ny"))
         {
             setNY(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-nz"))
         {
             setNZ(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-x"))
         {
             setXMinStat(atof(argv[i + 1]));
             setXMaxStat(atof(argv[i + 2]));
             i++; //requires two arguments
         }
         else if (!strcmp(argv[i], "-y"))
         {
             setYMinStat(atof(argv[i + 1]));
             setYMaxStat(atof(argv[i + 2]));
             i++; //requires two arguments
         }
         else if (!strcmp(argv[i], "-z"))
         {
             setZMinStat(atof(argv[i + 1]));
             setZMaxStat(atof(argv[i + 2]));
             i++; //requires two arguments
         }
         else if (!strcmp(argv[i], "-positions"))
         {
             loadPositions(argv[i+1]);
         }
         else if (!strcmp(argv[i], "-indSpecies"))
         {
             indSpecies = atoi(argv[i + 1]);
         }
         else if (!strcmp(argv[i], "-rmin"))
         {
             rmin = atof(argv[i + 1]);
         }
         else if (!strcmp(argv[i], "-rmax"))
         {
             rmax = atof(argv[i + 1]);
         }
         else if (!strcmp(argv[i], "-hmax"))
         {
             hmax = atof(argv[i + 1]);
         }
         else if (!strcmp(argv[i], "-periodicwalls"))
         {
             setDoPeriodicWalls(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-walls"))
         {
             setCGWidthalls(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-verbosity"))
         {
             setVerbosityLevel(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-verbose"))
         {
             verbose();
             i--; //requires no argument
         }
         else if (!strcmp(argv[i], "-format"))
         {
             format = atoi(argv[i + 1]);
         }
         else if (!strcmp(argv[i], "-doDoublePoints"))
         {
             if (argc <= i + 1 || argv[i + 1][0] == '-')
             {
                 setDoDoublePoints(true);
                 i--; //no argument
             }
             else
                 setDoDoublePoints(atoi(argv[i + 1]));
             std::cout << "set doDoublePoints=" << getDoDoublePoints() << std::endl;
         }
         else if (!strcmp(argv[i], "-w_over_rmax"))
         {
             setCGWidth_over_rmax(atof(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-tmin"))
         {
             setCGTimeMin(atof(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-tmax"))
         {
             setTimeMaxStat(atof(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-tint"))
         {
             setCGTimeAveragingInterval(atof(argv[i + 1]));
             //~ } else if (!strcmp(argv[i],"-timesteps")) {
             //~ setCGTimeAveragingInterval(atof(argv[i+1])*getTimeStep());
         }
         else if (!strcmp(argv[i], "-stepsize"))
         {
             int stepSize = atoi(argv[i + 1]);
             setStepSize(stepSize);
         }
         else if (!strcmp(argv[i], "-loaddatafile"))
         {
             readDataFile(argv[i + 1]);
             setCGTimeMin(getTime());
             setTimeMaxStat(getTimeMax());
         }
         else if (!strcmp(argv[i], "-statfilename") | !strcmp(argv[i], "-o"))
         {
             statFile.setName(argv[i + 1]);
             // setName(argv[i+1]);
         }
         else if (!strcmp(argv[i], "-timeaverage"))
         {
             setDoTimeAverage(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-timeaveragereset"))
         {
             nTimeAverageReset = atoi(argv[i + 1]);
         }
         else if (!strcmp(argv[i], "-gradient"))
         {
             // use default argument if no argument is given
             if (argc <= i + 1 || argv[i + 1][0] == '-')
             {
                 setDoGradient(true);
                 i--; //no argument
             }
             else
                 setDoGradient(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-timevariance"))
         {
             // use default argument if no argument is given
             if (argc <= i + 1 || argv[i + 1][0] == '-')
             {
                 setDoVariance(true);
                 i--; //no argument
             }
             else
                 setDoVariance(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-superexact"))
         {
             // use default argument if no argument is given
             if (argc <= i + 1 || argv[i + 1][0] == '-')
             {
                 setSuperExact(true);
                 i--; //no argument
             }
             else
                 setSuperExact(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-ignoreFixedParticles"))
         {
             // use default argument if no argument is given
             if (argc <= i + 1 || argv[i + 1][0] == '-')
             {
                 setDoIgnoreFixedParticles(true);
                 i--; //no argument
             }
             else
                 setDoIgnoreFixedParticles(atoi(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-StressTypeForFixedParticles"))
         {
             StressTypeForFixedParticles = atoi(argv[i + 1]);
         }
         else if (!strcmp(argv[i], "-mirrorAtDomainBoundary"))
         {
             setMirrorAtDomainBoundary(atof(argv[i + 1]));
         }
         else if (!strcmp(argv[i], "-stattype") || !strcmp(argv[i], "-statType"))
         {
             //this was already used in Statistics()
         }
         else if (readNextArgument(i, argc, argv))
         {
             //~ std::cout << " (read as MD argument, not statistics)" << std::endl; 
         }
         else
         {
             std::cerr << "Error: option unknown: " << argv[i] << std::endl;
             exit(-1);
         }
     }
 }

template<StatType T>

void StatisticsVector< T >::reset_statistics ( )

virtual

Set all statistical variables to zero.

Definition at line 97 of file StatisticsVector.hcc.

 {
     for (unsigned int i = 0; i < Points.size(); i++)
         Points[i].set_zero();
     if (getDoGradient())
         for (unsigned int i = 0; i < Points.size(); i++)
         {
             dx[i].set_zero();
             dy[i].set_zero();
             dz[i].set_zero();
         }
     for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
     {
         (*it)->setPreviousPosition((*it)->getPosition());
     }
 }

template<StatType T>

bool StatisticsVector< T >::satisfiesInclusionCriteria ( BaseParticle * P )

protected

Definition at line 1976 of file StatisticsVector.hcc.

References BaseInteractable::getIndSpecies(), BaseInteractable::getPosition(), BaseParticle::getRadius(), and Vec3D::Z.

 {
     return P->getRadius() > rmin
             && P->getRadius() < rmax
             && P->getPosition().Z < hmax
             && (indSpecies<0 || P->getIndSpecies() == indSpecies);
 }

template<StatType T>

void StatisticsVector< T >::set_h ( Mdouble h )

inline

Definition at line 184 of file StatisticsVector.h.

Referenced by main().

     {
         set_hx(h);
         set_hy(h);
         set_hz(h);
     }

template<StatType T>

void StatisticsVector< T >::set_hmax ( Mdouble new_ )

inline

Definition at line 884 of file StatisticsVector.h.

     {
         hmax = new_;
     }

template<StatType T>

void StatisticsVector< T >::set_hx ( Mdouble hx )

inline

Definition at line 132 of file StatisticsVector.h.

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

     {
         setNX(static_cast<int>(std::ceil((getXMaxStat() - getXMinStat()) / hx)));
     }

template<StatType T>

void StatisticsVector< T >::set_hy ( Mdouble hy )

inline

Definition at line 148 of file StatisticsVector.h.

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

     {
         setNY(static_cast<int>(std::ceil((getYMaxStat() - getYMinStat()) / hy)));
     }

template<StatType T>

void StatisticsVector< T >::set_hz ( Mdouble hz )

inline

Definition at line 166 of file StatisticsVector.h.

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

     {
         setNZ(static_cast<int>(std::ceil((getZMaxStat() - getZMinStat()) / hz)));
     }

template<StatType T>

void StatisticsVector< T >::set_infiniteStressForFixedParticles ( bool new_ )

inline

Definition at line 427 of file StatisticsVector.h.

     {
         StressTypeForFixedParticles = 2 + new_;
     }

template<StatType T>

void StatisticsVector< T >::set_Polynomial	(	std::vector< Mdouble >	new_coefficients,
		unsigned int	new_dim
	)

inline

Definition at line 804 of file StatisticsVector.h.

     {
         CGPolynomial.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

Definition at line 812 of file StatisticsVector.h.

     {
         CGPolynomial.set_polynomial(new_coefficients, num_coeff, new_dim);
     }

template<StatType T>

void StatisticsVector< T >::set_rmax ( Mdouble new_ )

inline

Definition at line 876 of file StatisticsVector.h.

     {
         rmax = new_;
     }

template<StatType T>

void StatisticsVector< T >::set_rmin ( Mdouble new_ )

inline

Definition at line 868 of file StatisticsVector.h.

     {
         rmin = new_;
     }

template<StatType T>

void StatisticsVector< T >::set_statType ( )

template<>

void StatisticsVector< RAZ >::set_statType ( )

Definition at line 2539 of file StatisticsVector.hcc.

References RAZ.

 {
     statType = RAZ;
 }

template<>

void StatisticsVector< RZ >::set_statType ( )

Definition at line 2543 of file StatisticsVector.hcc.

References RZ.

 {
     statType = RZ;
 }

template<>

void StatisticsVector< AZ >::set_statType ( )

Definition at line 2547 of file StatisticsVector.hcc.

References AZ.

 {
     statType = AZ;
 }

template<>

void StatisticsVector< RA >::set_statType ( )

Definition at line 2551 of file StatisticsVector.hcc.

References RA.

 {
     statType = RA;
 }

template<>

void StatisticsVector< A >::set_statType ( )

Definition at line 2555 of file StatisticsVector.hcc.

References A.

 {
     statType = A;
 }

template<>

void StatisticsVector< R >::set_statType ( )

Definition at line 2559 of file StatisticsVector.hcc.

References R.

 {
     statType = R;
 }

template<>

void StatisticsVector< XYZ >::set_statType ( )

Definition at line 2563 of file StatisticsVector.hcc.

References XYZ.

 {
     statType = XYZ;
 }

template<>

void StatisticsVector< XY >::set_statType ( )

Definition at line 2567 of file StatisticsVector.hcc.

References XY.

 {
     statType = XY;
 }

template<>

void StatisticsVector< XZ >::set_statType ( )

Definition at line 2571 of file StatisticsVector.hcc.

References XZ.

 {
     statType = XZ;
 }

template<>

void StatisticsVector< YZ >::set_statType ( )

Definition at line 2575 of file StatisticsVector.hcc.

References YZ.

 {
     statType = YZ;
 }

template<>

void StatisticsVector< X >::set_statType ( )

Definition at line 2579 of file StatisticsVector.hcc.

References X.

 {
     statType = X;
 }

template<>

void StatisticsVector< Y >::set_statType ( )

Definition at line 2583 of file StatisticsVector.hcc.

References Y.

 {
     statType = Y;
 }

template<>

void StatisticsVector< Z >::set_statType ( )

Definition at line 2587 of file StatisticsVector.hcc.

References Z.

 {
     statType = Z;
 }

template<StatType T>

void StatisticsVector< T >::setCGShape ( const char * new_ )

Definition at line 726 of file StatisticsVector.hcc.

References Gaussian, HeavisideSphere, and Polynomial.

Referenced by main().

 {
     if (!strcmp(new_, "Heaviside"))
     {
         int dim = 3;
         int numcoeff = 1;
         double heavisidecoeff[] = {1};
         set_Polynomial(heavisidecoeff, numcoeff, dim);
     }
     else if (!strcmp(new_, "Linear"))
     {
         int dim = 3;
         int numcoeff = 2;
         double lincoeff[] = {-1, 1};
         set_Polynomial(lincoeff, numcoeff, dim);
     }
     else if (!strcmp(new_, "Lucy"))
     {
         int dim = 3;
         int numcoeff = 5;
         double lucycoeff[] = {-3, 8, -6, 0, 1};
         set_Polynomial(lucycoeff, numcoeff, dim);
     }
     else if (!strcmp(new_, "Gaussian"))
     {
         setCGShape(Gaussian);
         return;
     }
     else if (!strcmp(new_, "HeavisideSphere"))
     {
         setCGShape(HeavisideSphere);
         return;
     }
     else
     {
         std::cerr << "Error: unknown CG_type: " << new_ << std::endl;
         exit(-1);
     }
     setPolynomialName(new_);
     setCGShape(Polynomial);
 }

template<StatType T>

void StatisticsVector< T >::setCGShape ( CG new_ )

Definition at line 769 of file StatisticsVector.hcc.

References Polynomial.

 {
     if (new_ == Polynomial && CGPolynomial.getOrder() < 0)
     {
         std::cerr << "Polynomial is not specified" << std::endl;
         exit(-1);
     }
     CG_type = new_;
 }

template<StatType T>

void StatisticsVector< T >::setCGTimeAveragingInterval ( Mdouble t )

inline

Definition at line 234 of file StatisticsVector.h.

     {
         tintStat = t;
     }

template<StatType T>

void StatisticsVector< T >::setCGTimeMin ( Mdouble t )

inline

Definition at line 218 of file StatisticsVector.h.

Referenced by main().

     {
         tminStat = t;
     }

template<StatType T>

void StatisticsVector< T >::setCGWidth ( Mdouble w )

inline

Set CG variables w2 and CG_invvolume.

Definition at line 312 of file StatisticsVector.h.

Referenced by main().

     {
         setCGWidth2(mathsFunc::square(w));
     }

template<StatType T>

void StatisticsVector< T >::setCGWidth2 ( Mdouble new_ )

Set CG variables w2 and CG_invvolume.

Todo:: {a default cutoff of 4.0 would be nicer}

Definition at line 67 of file StatisticsVector.hcc.

References Gaussian, HeavisideSphere, Polynomial, and mathsFunc::square().

Referenced by StatisticsVector< T >::initialiseStatistics(), and StatisticsVector< O >::setCGWidth().

 {
     // set w2
     if (new_ > 0.0)
         w2 = new_;
     else
         w2 = mathsFunc::square(w_over_rmax * particleHandler.getLargestParticle()->getRadius());
     //(2.0*getXMax()/nx);
     
     // set cutoff radius
     if (getCGShape() == HeavisideSphere || getCGShape() == Polynomial)
     {
         cutoff = sqrt(w2);
     }
     else if (getCGShape() == Gaussian)
     {
         if (getSuperExact())
             cutoff = 5.0 * sqrt(w2);
         else
             cutoff = 3.0 * sqrt(w2);
     }
     else
     {
         std::cerr << "error in CG_function" << std::endl;
         exit(-1);
     }
 }

template<StatType T>

void StatisticsVector< T >::setCGWidth_over_rmax ( Mdouble new_ )

inline

Definition at line 575 of file StatisticsVector.h.

     {
         w_over_rmax = new_;
     }

template<StatType T>

void StatisticsVector< T >::setCGWidthalls ( bool new_ )

inline

Definition at line 542 of file StatisticsVector.h.

     {
         walls = new_;
     }

template<StatType T>

void StatisticsVector< T >::setDoDoublePoints ( bool new_ )

inline

Definition at line 836 of file StatisticsVector.h.

     {
         doDoublePoints = new_;
     }

template<StatType T>

void StatisticsVector< T >::setDoGradient ( bool new_ )

inline

Definition at line 470 of file StatisticsVector.h.

     {
         doGradient = new_;
     }

template<StatType T>

void StatisticsVector< T >::setDoIgnoreFixedParticles ( bool new_ )

inline

Definition at line 502 of file StatisticsVector.h.

     {
         ignoreFixedParticles = new_;
     }

template<StatType T>

void StatisticsVector< T >::setDoPeriodicWalls ( bool new_ )

inline

Definition at line 558 of file StatisticsVector.h.

Referenced by main(), and Statistics().

     {
         periodicWalls = new_;
     }

template<StatType T>

void StatisticsVector< T >::setDoTimeAverage ( bool new_ )

inline

Definition at line 392 of file StatisticsVector.h.

     {
         doTimeAverage = new_;
     }

template<StatType T>

void StatisticsVector< T >::setDoVariance ( bool new_ )

inline

Definition at line 453 of file StatisticsVector.h.

     {
         doVariance = new_;
     }

template<StatType T>

Mdouble StatisticsVector< T >::setInfinitelyLongDistance ( )

template<>

Mdouble StatisticsVector< RAZ >::setInfinitelyLongDistance ( )

todo{check}

Definition at line 2596 of file StatisticsVector.hcc.

 {
     return std::min(std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
             std::max((getXMaxStat() - P2.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P2.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X)),
             std::max((getYMaxStat() - P2.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P2.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y));
 }

template<>

Mdouble StatisticsVector< RA >::setInfinitelyLongDistance ( )

Definition at line 2603 of file StatisticsVector.hcc.

 {
     return std::min(std::max((getYMaxStat() - P2.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P2.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y),
             std::max((getXMaxStat() - P2.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P2.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X));
 }

template<>

Mdouble StatisticsVector< RZ >::setInfinitelyLongDistance ( )

Definition at line 2608 of file StatisticsVector.hcc.

 {
     return std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
             std::max((getXMaxStat() - P2.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P2.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X));
 }

template<>

Mdouble StatisticsVector< AZ >::setInfinitelyLongDistance ( )

Definition at line 2613 of file StatisticsVector.hcc.

 {
     return std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
             std::max((getYMaxStat() - P2.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P2.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y));
 }

template<>

Mdouble StatisticsVector< R >::setInfinitelyLongDistance ( )

Definition at line 2618 of file StatisticsVector.hcc.

 {
     return std::max((getXMaxStat() - P2.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P2.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X);
 }

template<>

Mdouble StatisticsVector< A >::setInfinitelyLongDistance ( )

Definition at line 2622 of file StatisticsVector.hcc.

 {
     return std::max((getYMaxStat() - P2.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P2.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y);
 }

template<>

Mdouble StatisticsVector< XYZ >::setInfinitelyLongDistance ( )

Definition at line 2627 of file StatisticsVector.hcc.

 {
     return std::min(std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
             std::max((getXMaxStat() - P1.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P1.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X)),
             std::max((getYMaxStat() - P1.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P1.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y));
 }

template<>

Mdouble StatisticsVector< XY >::setInfinitelyLongDistance ( )

Definition at line 2633 of file StatisticsVector.hcc.

 {
     return std::min(std::max((getYMaxStat() - P1.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P1.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y),
             std::max((getXMaxStat() - P1.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P1.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X));
 }

template<>

Mdouble StatisticsVector< XZ >::setInfinitelyLongDistance ( )

Definition at line 2638 of file StatisticsVector.hcc.

 {
     return std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
             std::max((getXMaxStat() - P1.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P1.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X));
 }

template<>

Mdouble StatisticsVector< YZ >::setInfinitelyLongDistance ( )

Definition at line 2643 of file StatisticsVector.hcc.

 {
     return std::min(fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z),
             std::max((getYMaxStat() - P1.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P1.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y));
 }

template<>

Mdouble StatisticsVector< X >::setInfinitelyLongDistance ( )

Definition at line 2648 of file StatisticsVector.hcc.

 {
     return std::max((getXMaxStat() - P1.X + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X, -(P1.X - getXMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.X);
 }

template<>

Mdouble StatisticsVector< Y >::setInfinitelyLongDistance ( )

Definition at line 2652 of file StatisticsVector.hcc.

 {
     return std::max((getYMaxStat() - P1.Y + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y, -(P1.Y - getYMinStat() + 5 * sqrt(getCGWidthSquared())) / P1_P2_normal.Y);
 }

template<>

Mdouble StatisticsVector< Z >::setInfinitelyLongDistance ( )

Definition at line 2656 of file StatisticsVector.hcc.

 {
     return fabs((P1.Z - getZMinStat() + 5 * getCGWidth()) / P1_P2_normal.Z);
 }

template<StatType T>

void StatisticsVector< T >::setMirrorAtDomainBoundary ( Mdouble new_ )

inline

Definition at line 436 of file StatisticsVector.h.

     {
         mirrorAtDomainBoundary = new_;
     }

template<StatType T>

void StatisticsVector< T >::setN ( int n )

inline

Definition at line 174 of file StatisticsVector.h.

Referenced by main().

     {
         setNX(n);
         setNY(n);
         setNZ(n);
     }

template<StatType T>

void StatisticsVector< T >::setN	(	int	nx_,
		int	ny_,
		int	nz_
	)

inline

Definition at line 296 of file StatisticsVector.h.

     {
         nx = nx_;
         ny = ny_;
         nz = nz_;
     }

template<StatType T>

void StatisticsVector< T >::setNTimeAverageReset ( int new_ )

inline

Definition at line 852 of file StatisticsVector.h.

     {
         nTimeAverageReset = new_;
     }

template<StatType T>

void StatisticsVector< T >::setNX ( int new_ )

inline

Definition at line 124 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::set_hx(), and StatisticsVector< O >::setN().

     {
         nx = new_;
     }

template<>

void StatisticsVector< YZ >::setNX ( int new_ UNUSED )

Definition at line 2498 of file StatisticsVector.hcc.

 {
     nx = 1;
 }

template<>

void StatisticsVector< Y >::setNX ( int new_ UNUSED )

Definition at line 2510 of file StatisticsVector.hcc.

 {
     nx = 1;
 }

template<>

void StatisticsVector< Z >::setNX ( int new_ UNUSED )

Definition at line 2518 of file StatisticsVector.hcc.

 {
     nx = 1;
 }

template<StatType T>

void StatisticsVector< T >::setNY ( int new_ )

inline

Definition at line 140 of file StatisticsVector.h.

Referenced by main(), StatisticsVector< O >::set_hy(), and StatisticsVector< O >::setN().

     {
         ny = new_;
     }

template<>

void StatisticsVector< XZ >::setNY ( int new_ UNUSED )

Definition at line 2494 of file StatisticsVector.hcc.

 {
     ny = 1;
 }

template<>

void StatisticsVector< X >::setNY ( int new_ UNUSED )

Definition at line 2502 of file StatisticsVector.hcc.

 {
     ny = 1;
 }

template<>

void StatisticsVector< Z >::setNY ( int new_ UNUSED )

Definition at line 2522 of file StatisticsVector.hcc.

 {
     ny = 1;
 }

template<StatType T>

void StatisticsVector< T >::setNZ ( int new_ )

inline

Definition at line 156 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::set_hz(), and StatisticsVector< O >::setN().

     {
         nz = new_;
         if (getSystemDimensions() < 3)
             std::cerr << "Warning in set_nz: dimension less than 3" << std::endl;
     }

template<>

void StatisticsVector< XY >::setNZ ( int new_ UNUSED )

Definition at line 2490 of file StatisticsVector.hcc.

 {
     nz = 1;
 }

template<>

void StatisticsVector< X >::setNZ ( int new_ UNUSED )

Definition at line 2506 of file StatisticsVector.hcc.

 {
     nz = 1;
 }

template<>

void StatisticsVector< Y >::setNZ ( int new_ UNUSED )

Definition at line 2514 of file StatisticsVector.hcc.

 {
     nz = 1;
 }

template<StatType T>

void StatisticsVector< T >::setPolynomialName ( const char * new_name )

inline

Definition at line 820 of file StatisticsVector.h.

     {
         CGPolynomial.setName(new_name);
     }

template<StatType T>

void StatisticsVector< T >::setPositions ( )

Set position of StatisticsPoint points and set variables to 0.

Set position of statistics points and set variables to 0.

Definition at line 1763 of file StatisticsVector.hcc.

References A, AZ, R, RA, RAZ, RZ, X, Vec3D::X, XY, XYZ, XZ, Y, Vec3D::Y, YZ, Z, and Vec3D::Z.

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

 {
     std::cout << "set " << " positions" << std::endl;
     int N = positions_.size();
     if (N!=0)
     {
         std::cout << "setting " << N << " positions" << std::endl;
         Points.resize(N);
         for (int i = 0; i<N; i++) {
             Points[i].setPosition(positions_[i]);
         }
         nx = N; ny=1; nz=1;
     }
     else //set automatically on mesh
     {
         //Dimensions of the domain / n
         Vec3D diff = Vec3D(
             (getXMaxStat() - getXMinStat()) / nx,
             (getYMaxStat() - getYMinStat()) / ny,
             (getZMaxStat() - getZMinStat()) / nz);
         //Center of the domain
         Vec3D avg = 0.5 * Vec3D(
             getXMaxStat() + getXMinStat(),
             getYMaxStat() + getYMinStat(),
             getZMaxStat() + getZMinStat());
         //Left endpoint of the domain
         Vec3D Min = Vec3D(
             getXMinStat(),
             getYMinStat(),
             getZMinStat());
 
         int num[3] = {0, 0, 0};
         if (getMirrorAtDomainBoundary() != 0.0)
         {
             if (statType == X || statType == XY || statType == XZ || statType == XYZ)
             {
                 //add so many points to domain on both ends
                 //todo: the 8 is arbitrary
                 num[0] = static_cast<int>(std::floor(std::max(getMirrorAtDomainBoundary(), getCutoff()) / diff.X));
                 nx += 2 * num[0];
                 Min.X -= num[0] * diff.X;
             }
             if (statType == Y || statType == XY || statType == YZ || statType == XYZ)
             {
                 num[1] = static_cast<int>(std::floor(std::max(getMirrorAtDomainBoundary(), getCutoff()) / diff.Y));
                 ny += 2 * num[1];
                 Min.Y -= num[1] * diff.Y;
             }
             if (statType == Z || statType == XZ || statType == XZ || statType == XYZ)
             {
                 num[2] = static_cast<int>(std::floor(std::max(getMirrorAtDomainBoundary(), getCutoff()) / diff.Z));
                 nz += 2 * num[2];
                 Min.Z -= num[2] * diff.Z;
             }
         }
 
         //Set position of statistics points and set variables to 0
         N = std::max(nx, 1) * std::max(ny, 1) * std::max(nz, 1);
         Points.resize(N);
         {
             int n = 0;
             for (int i = 0; i < std::max(nx, 1); i++)
                 for (int j = 0; j < std::max(ny, 1); j++)
                     for (int k = 0; k < std::max(nz, 1); k++)
                     {
                         Points[n].setPosition(Vec3D((nx > 1) ? (Min.X + diff.X * (i + 0.5)) : avg.X,
                             (ny > 1) ? (Min.Y + diff.Y * (j + 0.5)) : avg.Y,
                             (nz > 1) ? (Min.Z + diff.Z * (k + 0.5)) : avg.Z));
                         if (doDoublePoints)
                             Points[n].setPosition(Points[n].getPosition()
                                 - Vec3D((i % 2) ? (0.99 * diff.X) : 0, (j % 2) ? (0.99 * diff.Y) : 0, (k % 2) ? (0.99 * diff.Z) : 0));
                         if (getMirrorAtDomainBoundary() != 0.0) if ((i < num[0]) || (i > nx - num[0]) || (j < num[1]) || (j > ny - num[1]) || (k < num[2]) || (k > nz - num[2]))
                         {
                             Points[n].mirrorParticle = n +
                                 +((i < num[0]) ? (2 * (num[0] - i) - 1) * ny * nz : 0) + ((i > nx - num[0]) ? (2 * (nx - num[0] - i) + 1) * ny * nz : 0)
                                 + ((j < num[1]) ? (2 * (num[1] - j) - 1) * nz : 0) + ((j > ny - num[1]) ? (2 * (ny - num[1] - j) + 1) * nz : 0)
                                 + ((k < num[2]) ? (2 * (num[2] - k) - 1) : 0) + ((k > nz - num[2]) ? (2 * (nz - num[2] - k) + 1) : 0);
                         }
                         n++;
                     }
         }
         if (statType == RAZ || statType == RZ || statType == RA || statType == AZ || statType == R || statType == A)
         {
             for (unsigned int k = 0; k < Points.size(); k++)
             {
                 Points[k].setPosition(Points[k].getPosition().getFromCylindricalCoordinates());
             }
         }
     }
     for (int i = 0; i<N; i++)
     {
         Points[i].setCGInverseVolume();
         Points[i].set_zero();
     }
     if (getDoGradient())
     {
         dx.resize(N);
         dy.resize(N);
         dz.resize(N);
         for (int n = 0; n < N; n++)
         {
             dx[n].setPosition(Points[n].getPosition());
             dx[n].mirrorParticle = Points[n].mirrorParticle;
             dx[n].set_zero();
             dy[n].setPosition(Points[n].getPosition());
             dy[n].mirrorParticle = Points[n].mirrorParticle;
             dy[n].set_zero();
             dz[n].setPosition(Points[n].getPosition());
             dz[n].mirrorParticle = Points[n].mirrorParticle;
             dz[n].set_zero();
         }
     }
 
     //Set TimeAverage the same way
     if (getDoTimeAverage())
     {
         timeAverage.resize(N);
         {
             int n = 0;
             for (int i = 0; i < std::max(nx, 1); i++)
                 for (int j = 0; j < std::max(ny, 1); j++)
                     for (int k = 0; k < std::max(nz, 1); k++)
                     {
                         timeAverage[n].setPosition(Points[n].getPosition());
                         timeAverage[n].set_zero();
                         timeAverage[n].mirrorParticle = Points[n].mirrorParticle;
                         n++;
                     }
         }
 
         if (getDoVariance())
         {
             timeVariance.resize(N);
             for (int n = 0; n < N; n++)
             {
                 timeVariance[n].setPosition(timeAverage[n].getPosition());
                 timeVariance[n].mirrorParticle = Points[n].mirrorParticle;
                 timeVariance[n].set_zero();
             }
         }
 
         if (getDoGradient())
         {
             dxTimeAverage.resize(N);
             dyTimeAverage.resize(N);
             dzTimeAverage.resize(N);
             for (int n = 0; n < N; n++)
             {
                 dxTimeAverage[n].setPosition(timeAverage[n].getPosition());
                 dxTimeAverage[n].mirrorParticle = Points[n].mirrorParticle;
                 dxTimeAverage[n].set_zero();
                 dyTimeAverage[n].setPosition(timeAverage[n].getPosition());
                 dyTimeAverage[n].mirrorParticle = Points[n].mirrorParticle;
                 dyTimeAverage[n].set_zero();
                 dzTimeAverage[n].setPosition(timeAverage[n].getPosition());
                 dzTimeAverage[n].mirrorParticle = Points[n].mirrorParticle;
                 dzTimeAverage[n].set_zero();
             }
         }
 
     } //end if (getDoTimeAverage())
 }

template<StatType T>

void StatisticsVector< T >::setStepSize ( unsigned int stepSize )

inline

Definition at line 924 of file StatisticsVector.h.

     {
         stepSize_ = stepSize;
     }

template<StatType T>

void StatisticsVector< T >::setStressTypeForFixedParticles ( int new_ )

inline

Definition at line 409 of file StatisticsVector.h.

     {
         StressTypeForFixedParticles = new_;
     }

template<StatType T>

void StatisticsVector< T >::setSuperExact ( bool new_ )

inline

Definition at line 486 of file StatisticsVector.h.

Referenced by main().

     {
         superexact = new_;
     }

template<StatType T>

void StatisticsVector< T >::setTimeMaxStat ( Mdouble t )

inline

Definition at line 226 of file StatisticsVector.h.

Referenced by main().

     {
         tmaxStat = t;
     }

template<StatType T>

void StatisticsVector< T >::setVerbosityLevel ( int new_ )

inline

Definition at line 526 of file StatisticsVector.h.

Referenced by main().

     {
         verbosity = new_;
     }

template<StatType T>

void StatisticsVector< T >::setXMaxStat ( Mdouble xmaxStat_ )

inline

Definition at line 770 of file StatisticsVector.h.

Referenced by main().

     {
         xmaxStat = xmaxStat_;
     }

template<StatType T>

void StatisticsVector< T >::setXMinStat ( Mdouble xminStat_ )

inline

Definition at line 746 of file StatisticsVector.h.

Referenced by main().

     {
         xminStat = xminStat_;
     }

template<StatType T>

void StatisticsVector< T >::setYMaxStat ( Mdouble ymaxStat_ )

inline

Definition at line 778 of file StatisticsVector.h.

Referenced by main().

     {
         ymaxStat = ymaxStat_;
     }

template<StatType T>

void StatisticsVector< T >::setYMinStat ( Mdouble yminStat_ )

inline

Definition at line 754 of file StatisticsVector.h.

Referenced by main().

     {
         yminStat = yminStat_;
     }

template<StatType T>

void StatisticsVector< T >::setZMaxStat ( Mdouble zmaxStat_ )

inline

Definition at line 786 of file StatisticsVector.h.

Referenced by main().

     {
         zmaxStat = zmaxStat_;
     }

template<StatType T>

void StatisticsVector< T >::setZMinStat ( Mdouble zminStat_ )

inline

Definition at line 762 of file StatisticsVector.h.

Referenced by main().

     {
         zminStat = zminStat_;
     }

template<StatType T>

void StatisticsVector< T >::statistics_from_fstat_and_data ( )

get StatisticsPoint

Todo:: : should be heaviest particle

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

Definition at line 1113 of file StatisticsVector.hcc.

References MULTIPLE_FILES, MULTIPLE_FILES_PADDED, ONE_FILE, and DPMBase::write().

Referenced by main(), and Statistics().

 {
     std::cout << "set " << " positions" << std::endl;
     dataFile.setCounter(0);
     fStatFile.setCounter(0);
 
     //This opens the file the data will be recalled from
     if (dataFile.getFileType() == FileType::ONE_FILE)
         dataFile.open(std::fstream::in);
     if (fStatFile.getFileType() == FileType::ONE_FILE)
         fStatFile.open(std::fstream::in);
 
 
     //load first set of timesteps to obtain xdim, ydim, zdim
     
     //try to read data files; if data.0000 does not exist, try the next, up to .9999
     for (int i = 1; i < 10000; i++)
     {
         if (readNextDataFile(format))
             break;
         fStatFile.openNextFile(std::fstream::in);
     }
 
     //set tminstat to smallest time evaluated
     if (getTime() >= getCGTimeMin())
         setCGTimeMin(getTime() - getTimeStep());
 
     //This creates the file statistics will be saved to
     //getStatFile.open(std::fstream::out);
     
     if (verbosity > 1)
     {
         std::cout << "Input from " << dataFile.getName() << std::endl;
         std::cout << "Input from " << fStatFile.getName() << std::endl;
         std::cout << "Output to " << statFile.getName() << std::endl << std::endl;
     }
     
     initialiseStatistics();
     
     //std::couts collision time
 //    if (verbosity > 1)
 //    {
 //        //getLargestParticle()->computeMass(getSpecies()); //This should already have been done in readNextdataFile;
 //        Mdouble mass = particleHandler.getLargestParticle()->getMass();
 //        std::cout << "Collision Time " << speciesHandler.getObject(0)->getCollisionTime(mass)
 //            << " and restitution coefficient " << speciesHandler.getObject(0)->getRestitutionCoefficient(mass)
 //                << " for mass " << mass << std::endl
 //                << std::endl;
 //    }
     //std::couts particles
     if (verbosity > 2)
     {
         DPMBase::write(std::cout, true);
         std::cout << std::endl;
     }
     else if (verbosity)
     {
         //DPMBase::write(std::cout, false);
         //std::cout << std::endl;
     }
 
 
     //this increases the file counter
     if ((fStatFile.getFileType() == FileType::MULTIPLE_FILES || fStatFile.getFileType() == FileType::MULTIPLE_FILES_PADDED) && getTime() < getCGTimeMin())
     {
         if (verbosity > 1)
             findNextExistingDataFile(getCGTimeMin());
         else
             findNextExistingDataFile(getCGTimeMin(), false);
         for (unsigned int i = 1; i < dataFile.getCounter(); i++)
             jump_fstat();
 
     }
 
     //Output statistics for each time step
     std::cout << "Start statistics" << std::endl;
     do
     {
         if (getTime() > getTimeMaxStat())
         {
             std::cout
                 << "reached end t=" << std::setprecision(9) << getTime()
                 << " tmaxStat=" << std::setprecision(9) << getTimeMaxStat()
                 << std::endl;
             break;
         }
         else if (getTime() < getCGTimeMin())
         {
             if (verbosity > 1)
                 std::cout << "Jumping statistics t=" << getTime() << " tminStat()=" << getCGTimeMin() << std::endl;
             jump_fstat();
         }
         else
         {
             if (verbosity)
                 std::cout << "Get statistics for t=" << std::setprecision(9) << getTime() << std::endl;
             if (verbosity > 1)
                 std::cout << "#particles=" << particleHandler.getNumberOfObjects() << std::endl;
             
             if (!walls)
                 wallHandler.clear();
             if (!periodicWalls)
                 boundaryHandler.clear();
             
             outputStatistics();
             gather_force_statistics_from_fstat_and_data();
             processStatistics(true);
             
         }
         //if step size>1, skip a number of steps.
         if (dataFile.getFileType() == FileType::ONE_FILE)
             for (unsigned int i = 1; i < getStepSize(); i++)
             {
                 //you only get here if you change the step size and you use a single data file
                 readNextDataFile(format);
                 jump_fstat();
             }
     } while (readNextDataFile(format));
     //set tmax to largest time evaluated
     if (getTime() < getTimeMaxStat())
         setTimeMaxStat(getTime());
     
     finishStatistics();
     dataFile.close();
     fStatFile.close();
 }

template<StatType T>

void StatisticsVector< T >::statistics_from_p3 ( )

this is a modified version of statistics_from_fstat_and_data. It is used to read p3d files

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

Definition at line 1343 of file StatisticsVector.hcc.

References helpers::openFile(), DPMBase::outputXBallsData(), DPMBase::setSystemDimensions(), DPMBase::write(), DPMBase::writeRestartFile(), and DPMBase::writeXBallsScript().

Referenced by main().

 {
     //This opens the files the data will be recalled from
     std::string p3p_filename = getName().c_str();
     p3p_filename += ".p3p";
     bool opened = helpers::openFile(p3p_file, p3p_filename, std::fstream::in);
     std::string p3c_filename = getName().c_str();
     p3c_filename += ".p3c";
     helpers::openFile(p3c_file, p3c_filename, std::fstream::in);
     std::string p3w_filename = getName().c_str();
     p3w_filename += ".p3w";
     helpers::openFile(p3w_file, p3w_filename, std::fstream::in);
     int version = 3;
     //If p3 is not available, open p4 files
     if (!opened)
     {
         std::cerr << "Warning could not open " << p3p_filename << std::endl;
         version = 4;
         p3p_filename = getName().c_str();
         p3p_filename += ".p4p";
         opened = helpers::openFile(p3p_file, p3p_filename, std::fstream::in);
         p3c_filename = getName().c_str();
         p3c_filename += ".p4c";
         opened = helpers::openFile(p3c_file, p3c_filename, std::fstream::in);
         p3w_filename = getName().c_str();
         p3w_filename += ".p4w";
         opened = helpers::openFile(p3w_file, p3w_filename, std::fstream::in);
         if (!opened)
         {
             std::cerr << "could not open " << p3p_filename << std::endl;
             exit(-1);
         }
     }
     std::cout << "opened " << p3p_filename << std::endl;
     
     //load first set of timesteps to obtain xdim, ydim, zdim
     
     //try to read data files; if data.0000 does not exist, try the next, up to .9999
     read_next_from_p3p_file();
     
     //set tminstat to smallest time evaluated
     if (getTime() >= getCGTimeMin())
         setCGTimeMin(getTime() - getTimeStep());
     DPMBase::writeRestartFile();
     
     dataFile.open(std::fstream::out);
     DPMBase::setSystemDimensions(3);
     setParticleDimensions(3);
     DPMBase::outputXBallsData(dataFile.getFstream());
     DPMBase::writeXBallsScript();
     //exit(-1);
     
     //This creates the file statistics will be saved to
     //getStatFile.open(std::fstream::out);
     
     initialiseStatistics();
 
     //std::couts collision time
     if (verbosity > 1)
     {
         //getLargestParticle()->computeMass(getSpecies()); //This should already have been done in read_next_from_p3p_file();
         Mdouble mass = particleHandler.getLargestParticle()->getMass();
         auto species = dynamic_cast<LinearViscoelasticSpecies*>(speciesHandler.getObject(0));
         if (species!= nullptr)
         {
             std::cout << "Collision Time " << species->getCollisionTime(mass)
                 << " and restitution coefficient " << species->getRestitutionCoefficient(mass)
                 << " for mass " << mass << std::endl
                 << std::endl;
         }
         //std::couts particles
     }
     if (verbosity > 2)
     {
         DPMBase::write(std::cout, true);
         std::cout << std::endl;
     }
     else if (verbosity)
     {
         DPMBase::write(std::cout,false);
         std::cout << std::endl;
     }
 
     //Output statistics for each time step
     std::cout << "Start statistics" << std::endl;
     do
     {
         if (getTime() > getTimeMaxStat())
         {
             std::cout << "reached end t=" << std::setprecision(9) << getTime() << " tmaxStat=" << std::setprecision(9) << getTimeMaxStat() << std::endl;
             break;
         }
         else if (getTime() >= getCGTimeMin())
         {
             if (verbosity)
                 std::cout << "Get statistics for t=" << std::setprecision(9) << getTime() << std::endl;
             if (verbosity > 1)
                 std::cout << "#particles=" << particleHandler.getNumberOfObjects() << std::endl;
             
             if (!walls)
                 wallHandler.clear();
             if (!periodicWalls)
                 boundaryHandler.clear();
             
             outputStatistics();
             gather_force_statistics_from_p3c(version);
             processStatistics(true);
             
         }
         else
         {
             if (verbosity > 1)
                 std::cout
                 << "Jumping statistics t=" << getTime()
                         << " tminStat()=" << getCGTimeMin() << std::endl;
             jump_p3c();
         }
         //if step size>1, skip a number of steps.
         for (unsigned int i = 1; i < getStepSize(); i++)
         {
             //you only get here if you change the step size and you use a single data file
             read_next_from_p3p_file();
             jump_p3c();
         }
     } while (read_next_from_p3p_file());
     //set tmax to largest time evaluated
     if (getTime() < getTimeMaxStat())
         setTimeMaxStat(getTime());
     
     finishStatistics();
     p3p_file.close();
     p3c_file.close();
     p3w_file.close();
 }

template<StatType T>

void StatisticsVector< T >::verbose ( )

inline

Definition at line 518 of file StatisticsVector.h.

     {
         verbosity = 2;
     }

template<StatType T>

void StatisticsVector< T >::write_statistics ( )

Writes regular statistics.

Definition at line 936 of file StatisticsVector.hcc.

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

template<StatType T>

void StatisticsVector< T >::write_time_average_statistics ( )

Writes out time averaged statistics.

Definition at line 1026 of file StatisticsVector.hcc.

 {
     static bool first = true;
     
     if (verbosity)
         std::cout << std::endl << "averaging " << nTimeAverage << " timesteps " << std::endl;
     if (nTimeAverage == 0)
     {
         fprintf(stderr, "\n\n\tERROR :: Can not do TimeAverage of 0 timesteps\n\n");
         //exit(-1);
         return;
     }
     if (first)
     {
         std::cout << "first" << nTimeAverage << std::endl;
         for (unsigned int i = 0; i < timeAverage.size(); i++)
         {
             timeAverage[i].firstTimeAverage(nTimeAverage);
             if (getDoVariance())
             {
                 timeVariance[i].firstTimeAverage(nTimeAverage);
                 timeVariance[i] -= timeAverage[i].getSquared();
             }
             if (getDoGradient())
             {
                 dxTimeAverage[i].firstTimeAverage(nTimeAverage);
                 dyTimeAverage[i].firstTimeAverage(nTimeAverage);
                 dzTimeAverage[i].firstTimeAverage(nTimeAverage);
             }
         }
         first = false;
     }
     else
     {
         std::cout << "next" << nTimeAverage << std::endl;
         for (unsigned int i = 0; i < timeAverage.size(); i++)
         {
             timeAverage[i] /= nTimeAverage;
             if (getDoVariance())
             {
                 timeVariance[i] /= nTimeAverage;
                 timeVariance[i] -= timeAverage[i].getSquared();
             }
             if (getDoGradient())
             {
                 dxTimeAverage[i] /= nTimeAverage;
                 dyTimeAverage[i] /= nTimeAverage;
                 dzTimeAverage[i] /= nTimeAverage;
             }
         }
     }
     
     // write average to .stat file
     statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getCGTimeMin() << " " << getTime() << std::endl;
     for (unsigned int i = 0; i < timeAverage.size(); i++)
         statFile.getFstream() << timeAverage[i];
     // write to std::cout
     StatisticsPoint<T> avg = average(timeAverage);
     std::cout << "Averages: " << avg.print() << std::endl;
     
     if (getDoVariance())
     {
         // write variance to .stat file
         statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getCGTimeMin() << " " << getTime() << " " << std::endl;
         for (unsigned int i = 0; i < timeVariance.size(); i++)
             statFile.getFstream() << timeVariance[i];
         //StatisticsPoint<T> var = average(timeVariance);
     } //end if (getDoVariance)
     
     if (getDoGradient())
     {
         // write variance to .stat file
         statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getCGTimeMin() << " " << getTime() << " " << std::endl;
         for (unsigned int i = 0; i < timeAverage.size(); i++)
             statFile.getFstream() << dxTimeAverage[i];
         statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getCGTimeMin() << " " << getTime() << " " << std::endl;
         for (unsigned int i = 0; i < timeAverage.size(); i++)
             statFile.getFstream() << dyTimeAverage[i];
         statFile.getFstream() << std::left << std::setprecision(8) << std::setw(6) << getCGTimeMin() << " " << getTime() << " " << std::endl;
         for (unsigned int i = 0; i < timeAverage.size(); i++)
             statFile.getFstream() << dzTimeAverage[i];
     }
     
     setCGTimeMin(getTime());
     
 }

template<StatType T>

void StatisticsVector< T >::writeOutputFiles ( )

virtual

Writes the simulation data onto all the files i.e. .data, .ene, .fstat ...

Reimplemented from DPMBase.

Definition at line 2666 of file StatisticsVector.hcc.

References DPMBase::gatherContactStatistics(), and DPMBase::writeOutputFiles().

 {
     if (statFile.saveCurrentTimestep(getNtimeSteps()))
     {
         outputStatistics();
         DPMBase::gatherContactStatistics();
         processStatistics(true);
     }
 
     DPMBase::writeOutputFiles();
 }

Member Data Documentation

template<StatType T>

CG StatisticsVector< T >::CG_type

protected

coarse graining type (Gaussian, Heaviside, Polynomial)

Definition at line 1025 of file StatisticsVector.h.

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

template<StatType T>

NORMALIZED_POLYNOMIAL<T> StatisticsVector< T >::CGPolynomial

protected

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

Definition at line 1029 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::evaluateIntegral(), StatisticsVector< O >::evaluatePolynomial(), StatisticsVector< O >::evaluatePolynomialGradient(), StatisticsVector< O >::getPolynomialName(), StatisticsVector< O >::set_Polynomial(), StatisticsVector< O >::setPolynomialName(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

Mdouble StatisticsVector< T >::cutoff

protected

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

Definition at line 1042 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getCutoff(), StatisticsVector< O >::getCutoff2(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

Mdouble StatisticsVector< T >::cutoff2

protected

Definition at line 1042 of file StatisticsVector.h.

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

template<StatType T>

bool StatisticsVector< T >::doDoublePoints

protected

Definition at line 1133 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getDoDoublePoints(), and StatisticsVector< O >::setDoDoublePoints().

template<StatType T>

bool StatisticsVector< T >::doGradient

protected

Determines if gradient is outputted.

Definition at line 1014 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getDoGradient(), StatisticsVector< O >::setDoGradient(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

bool StatisticsVector< T >::doTimeAverage

protected

Determines if output is averaged over time.

Definition at line 999 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getDoTimeAverage(), StatisticsVector< O >::setDoTimeAverage(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

bool StatisticsVector< T >::doVariance

protected

Determines if variance is outputted.

Definition at line 1009 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getDoVariance(), StatisticsVector< O >::setDoVariance(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

std::vector<StatisticsPoint<T> > StatisticsVector< T >::dx

protected

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

Definition at line 960 of file StatisticsVector.h.

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

template<StatType T>

std::vector<StatisticsPoint<T> > StatisticsVector< T >::dxTimeAverage

protected

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

Definition at line 984 of file StatisticsVector.h.

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

template<StatType T>

std::vector<StatisticsPoint<T> > StatisticsVector< T >::dy

protected

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

Definition at line 964 of file StatisticsVector.h.

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

template<StatType T>

std::vector<StatisticsPoint<T> > StatisticsVector< T >::dyTimeAverage

protected

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

Definition at line 989 of file StatisticsVector.h.

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

template<StatType T>

std::vector<StatisticsPoint<T> > StatisticsVector< T >::dz

protected

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

Definition at line 968 of file StatisticsVector.h.

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

template<StatType T>

std::vector<StatisticsPoint<T> > StatisticsVector< T >::dzTimeAverage

protected

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

Definition at line 994 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::format

protected

format of the data input file

Definition at line 1120 of file StatisticsVector.h.

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

template<StatType T>

Mdouble StatisticsVector< T >::hmax

protected

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

Definition at line 1085 of file StatisticsVector.h.

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

template<StatType T>

bool StatisticsVector< T >::ignoreFixedParticles

protected

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

Definition at line 1101 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getDoIgnoreFixedParticles(), StatisticsVector< O >::setDoIgnoreFixedParticles(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

Mdouble StatisticsVector< T >::indSpecies

protected

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

Definition at line 1069 of file StatisticsVector.h.

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

template<StatType T>

bool StatisticsVector< T >::isMDCLR

protected

Definition at line 1128 of file StatisticsVector.h.

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

template<StatType T>

Mdouble StatisticsVector< T >::mirrorAtDomainBoundary

protected

0: Statistics near the wall are equal to statistics away from the wall; 1: Statistics are mirrored at the domain boundaries; up to a maximum depth specified by this number

Definition at line 1125 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getMirrorAtDomainBoundary(), StatisticsVector< O >::setMirrorAtDomainBoundary(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

int StatisticsVector< T >::nTimeAverage

protected

A counter needed to average over time.

Definition at line 1019 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::nTimeAverageReset

protected

Determines after how many timesteps the time average is reset.

Definition at line 1004 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getNTimeAverageReset(), and StatisticsVector< O >::setNTimeAverageReset().

template<StatType T>

int StatisticsVector< T >::nx

protected

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].

Definition at line 940 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getN(), StatisticsVector< O >::getNX(), StatisticsVector< O >::setN(), StatisticsVector< O >::setNX(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

int StatisticsVector< T >::nxMirrored

protected

extension of grid size from mirrored points

Definition at line 950 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::ny

protected

see nx

Definition at line 941 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getN(), StatisticsVector< O >::getNY(), StatisticsVector< O >::setN(), StatisticsVector< O >::setNY(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

int StatisticsVector< T >::nyMirrored

protected

Definition at line 950 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::nz

protected

see nx

Definition at line 942 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getN(), StatisticsVector< O >::getNZ(), StatisticsVector< O >::setN(), StatisticsVector< O >::setNZ(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

int StatisticsVector< T >::nzMirrored

protected

Definition at line 950 of file StatisticsVector.h.

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

template<StatType T>

Vec3D StatisticsVector< T >::P1

protected

Position of first contact point.

Definition at line 1142 of file StatisticsVector.h.

template<StatType T>

Vec3D StatisticsVector< T >::P1_P2_CollisionalHeatFlux

protected

not yet working

Definition at line 1164 of file StatisticsVector.h.

template<StatType T>

Vec3D StatisticsVector< T >::P1_P2_Contact

protected

Definition at line 1154 of file StatisticsVector.h.

template<StatType T>

Matrix3D StatisticsVector< T >::P1_P2_ContactCoupleStress

protected

Definition at line 1152 of file StatisticsVector.h.

template<StatType T>

Mdouble StatisticsVector< T >::P1_P2_Dissipation

protected

not yet working

Definition at line 1166 of file StatisticsVector.h.

template<StatType T>

Mdouble StatisticsVector< T >::P1_P2_distance

protected

Length of contact line.

Definition at line 1148 of file StatisticsVector.h.

template<StatType T>

MatrixSymmetric3D StatisticsVector< T >::P1_P2_Fabric

protected

Fabric.

Definition at line 1162 of file StatisticsVector.h.

template<StatType T>

Vec3D StatisticsVector< T >::P1_P2_normal

protected

Direction of contact.

Definition at line 1146 of file StatisticsVector.h.

template<StatType T>

Matrix3D StatisticsVector< T >::P1_P2_NormalStress

protected

Contact stress from normal forces along the line of contact.

Definition at line 1150 of file StatisticsVector.h.

template<StatType T>

Vec3D StatisticsVector< T >::P1_P2_NormalTraction

protected

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

Definition at line 1158 of file StatisticsVector.h.

template<StatType T>

Mdouble StatisticsVector< T >::P1_P2_Potential

protected

not yet working

Definition at line 1168 of file StatisticsVector.h.

template<StatType T>

Matrix3D StatisticsVector< T >::P1_P2_TangentialStress

protected

Contact stress from tangential forces along the line of contact.

Definition at line 1156 of file StatisticsVector.h.

template<StatType T>

Vec3D StatisticsVector< T >::P1_P2_TangentialTraction

protected

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

Definition at line 1160 of file StatisticsVector.h.

template<StatType T>

Vec3D StatisticsVector< T >::P2

protected

Position of second contact point.

Definition at line 1144 of file StatisticsVector.h.

template<StatType T>

std::fstream StatisticsVector< T >::p3c_file

protected

Definition at line 1197 of file StatisticsVector.h.

template<StatType T>

std::fstream StatisticsVector< T >::p3p_file

protected

Definition at line 1196 of file StatisticsVector.h.

template<StatType T>

std::fstream StatisticsVector< T >::p3w_file

protected

Definition at line 1198 of file StatisticsVector.h.

template<StatType T>

bool StatisticsVector< T >::periodicWalls

protected

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

Todo:: {Thomas: the case periodicWalls=true seems to mess up some statistics. Needs to be checked or removed}

Definition at line 1096 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getDoPeriodicWalls(), StatisticsVector< O >::setDoPeriodicWalls(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

std::vector<StatisticsPoint<T> > StatisticsVector< T >::Points

protected

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

Definition at line 956 of file StatisticsVector.h.

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

template<StatType T>

std::vector<Vec3D> StatisticsVector< T >::positions_

protected

Definition at line 1212 of file StatisticsVector.h.

template<StatType T>

Mdouble StatisticsVector< T >::rmax

protected

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

Definition at line 1080 of file StatisticsVector.h.

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

template<StatType T>

Mdouble StatisticsVector< T >::rmin

protected

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()

Definition at line 1075 of file StatisticsVector.h.

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

template<StatType T>

StatType StatisticsVector< T >::statType

protected

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

Definition at line 936 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::set_statType().

template<StatType T>

unsigned int StatisticsVector< T >::stepSize_

protected

Definition at line 1210 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getStepSize(), and StatisticsVector< O >::setStepSize().

template<StatType T>

int StatisticsVector< T >::StressTypeForFixedParticles

protected

Stress type for fixed particles.

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

Definition at line 1110 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getStressTypeForFixedParticles(), StatisticsVector< O >::set_infiniteStressForFixedParticles(), StatisticsVector< O >::setStressTypeForFixedParticles(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

bool StatisticsVector< T >::superexact

protected

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

Definition at line 1130 of file StatisticsVector.h.

Referenced by StatisticsVector< O >::getSuperExact(), StatisticsVector< O >::setSuperExact(), and StatisticsVector< T >::StatisticsVector().

template<StatType T>

std::vector<StatisticsPoint<T> > StatisticsVector< T >::timeAverage

protected

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