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 ()
	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 (unsigned int argc, char *argv[])
	Advanced constructor that accepts arguments from the command line. More...

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

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

void	set_statType ()

void	print_help ()

void	set_nx (int new_)

void	set_hx (Mdouble hx)

void	set_ny (int new_)

void	set_hy (Mdouble hy)

void	set_nz (int new_)

void	set_hz (Mdouble hz)

void	set_n (int n)

void	set_h (Mdouble h)

int	get_nx ()

int	get_ny ()

int	get_nz ()

void	set_tminStat (Mdouble t)

void	set_tmaxStat (Mdouble t)

void	set_tintStat (Mdouble t)

Mdouble	get_tminStat ()

Mdouble	get_tmaxStat ()

Mdouble	get_tintStat ()

bool	check_current_time_for_statistics ()

void	set_CG_type (const char *new_)

void	set_CG_type (CG new_)

CG	get_CG_type ()

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

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

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

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

Mdouble	get_w ()

Mdouble	get_w2 ()

Mdouble	get_cutoff ()

Mdouble	get_cutoff2 ()

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	set_doTimeAverage (bool new_)

bool	get_doTimeAverage ()

void	set_StressTypeForFixedParticles (int new_)

int	get_StressTypeForFixedParticles ()

void	set_infiniteStressForFixedParticles (bool new_)

void	set_mirrorAtDomainBoundary (Mdouble new_)

Mdouble	get_mirrorAtDomainBoundary ()

void	set_doVariance (bool new_)

bool	get_doVariance ()

void	set_doGradient (bool new_)

bool	get_doGradient ()

void	set_superexact (bool new_)

bool	get_superexact ()

void	set_ignoreFixedParticles (bool new_)

bool	get_ignoreFixedParticles ()

void	verbose ()

void	set_verbosity (int new_)

int	get_verbosity ()

void	set_walls (bool new_)

bool	get_walls ()

void	set_periodicWalls (bool new_)

bool	get_periodicWalls ()

void	set_w_over_rmax (Mdouble new_)

Mdouble	get_w_over_rmax ()

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

bool	read_next_from_data_file (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	initialize_statistics ()
	Initializes statistics, i.e. setting w2, setting the grid and writing the header lines in the .stat file. More...

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

void	gather_statistics_collision (int index1, int index2, Vec3D Contact, Mdouble delta, Mdouble ctheta, Mdouble fdotn, Mdouble fdott, Vec3D P1_P2_normal_, Vec3D P1_P2_tangential)
	Calculates statistics for one collision (can be any kind of collision) More...

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

void	finish_statistics ()
	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 > >	get_Points ()

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

Mdouble	get_yminStat ()

Mdouble	get_zminStat ()

Mdouble	get_xmaxStat ()

Mdouble	get_ymaxStat ()

Mdouble	get_zmaxStat ()

void	set_xminStat (Mdouble xminStat_)

void	set_yminStat (Mdouble yminStat_)

void	set_zminStat (Mdouble zminStat_)

void	set_xmaxStat (Mdouble xmaxStat_)

void	set_ymaxStat (Mdouble ymaxStat_)

void	set_zmaxStat (Mdouble zmaxStat_)

int	get_nTimeAverage ()

Mdouble	setInfinitelyLongDistance ()

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

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

void	set_PolynomialName (const char *new_name)

std::string	get_PolynomialName ()

void	set_doublePoints (bool new_)

bool	get_doublePoints ()

void	set_TimeAverageReset (int new_)

bool	get_TimeAverageReset ()

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)

template<>
void	set_nz (int new_ UNUSED)

template<>
void	set_ny (int new_ UNUSED)

template<>
void	set_nx (int new_ UNUSED)

template<>
void	set_ny (int new_ UNUSED)

template<>
void	set_nz (int new_ UNUSED)

template<>
void	set_nx (int new_ UNUSED)

template<>
void	set_nz (int new_ UNUSED)

template<>
void	set_nx (int new_ UNUSED)

template<>
void	set_ny (int new_ UNUSED)

template<>
void	set_nx (int new_ UNUSED)

template<>
void	set_ny (int new_ UNUSED)

template<>
void	set_nz (int new_ UNUSED)

template<>
void	set_statType ()

template<>
void	set_statType ()

template<>
void	set_statType ()

template<>
void	set_statType ()

template<>
void	set_statType ()

template<>
void	set_statType ()

template<>
void	set_statType ()

template<>
void	set_statType ()

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

template<>
double	setInfinitelyLongDistance ()

Public Member Functions inherited from MD
void	constructor ()
	A public constructor, which sets defaults so the problem can be solved off the shelf. More...

	MD ()

	MD (STD_save &other)

virtual	~MD ()

void	info ()
	Set up a virtual info this will be provided from the inhertiance. More...

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

void	solve (unsigned int argc, char *argv[])
	Read arguments before solving. More...

void	solveWithMDCLR ()
	Tries to solve the problem using MDCLR. More...

Mdouble	get_t ()
	Access function for the time. More...

void	set_t (Mdouble new_t)
	Access function for the time. More...

int	get_NSpecies ()
	Allows the number of Species to be accessed. More...

std::vector< CSpecies > &	get_Species (void)
	Allows the species to be copied. More...

CSpecies *	get_Species (int i)
	Allows the species to be accessed. More...

CSpecies *	get_MixedSpecies (int i, int j)
	Allows the mixed species to be accessed. More...

void	set_MixedSpecies (int i, int j, CSpecies &S)
	Allows the mixed species to be set. More...

void	set_tmax (Mdouble new_tmax)
	Allows the upper time limit to be changed. More...

Mdouble	get_tmax ()
	Allows the upper time limit to be accessed. More...

ParticleHandler &	getParticleHandler ()

WallHandler &	getWallHandler ()

BoundaryHandler &	getBoundaryHandler ()

void	set_savecount (int new_)
	Allows the number of time steps between saves to be changed, see also set_number_of_saves. More...

void	set_save_count_all (int new_)

void	set_save_count_data (int new_)

void	set_save_count_ene (int new_)

void	set_save_count_stat (int new_)

void	set_save_count_fstat (int new_)

int	get_savecount ()
	Allows the number of time steps between saves to be accessed. More...

int	get_save_count ()

int	get_save_count_data ()

int	get_save_count_ene ()

int	get_save_count_stat ()

int	get_save_count_fstat ()

void	set_do_stat_always (bool new_)
	Sets how often the data is saved using the number of saves wanted, tmax, and dt. See also set_savecount. More...

void	set_number_of_saves (Mdouble N)

void	set_number_of_saves_all (Mdouble N)

void	set_number_of_saves_data (Mdouble N)

void	set_number_of_saves_ene (Mdouble N)

void	set_number_of_saves_stat (Mdouble N)

void	set_number_of_saves_fstat (Mdouble N)

void	set_plastic_k1_k2max_kc_depth (Mdouble k1_, Mdouble k2max_, Mdouble kc_, Mdouble depth_, unsigned int indSpecies=0)
	Allows the plastic constants to be changed. More...

void	set_k1 (Mdouble new_, unsigned int indSpecies=0)

void	set_k2max (Mdouble new_, unsigned int indSpecies=0)

void	set_kc (Mdouble new_, unsigned int indSpecies=0)

void	set_depth (Mdouble new_, unsigned int indSpecies=0)

Mdouble	get_k1 (unsigned int indSpecies=0)
	Allows the plastic constants to be accessed. More...

Mdouble	get_k2max (unsigned int indSpecies=0)

Mdouble	get_kc (unsigned int indSpecies=0)

Mdouble	get_depth (unsigned int indSpecies=0)

Mdouble	get_plastic_dt (Mdouble mass, unsigned int indSpecies=0)

void	set_k (Mdouble new_, unsigned int indSpecies=0)
	Allows the spring constant to be changed. More...

Mdouble	get_k (int indSpecies=0)
	Allows the spring constant to be accessed. More...

void	set_kt (Mdouble new_, unsigned int indSpecies=0)
	Allows the spring constant to be changed. More...

Mdouble	get_kt (int indSpecies=0)
	Allows the spring constant to be accessed. More...

void	set_krolling (Mdouble new_, unsigned int indSpecies=0)
	Allows the spring constant to be changed. More...

Mdouble	get_krolling (int indSpecies=0)
	Allows the spring constant to be accessed. More...

void	set_ktorsion (Mdouble new_, unsigned int indSpecies=0)
	Allows the spring constant to be changed. More...

Mdouble	get_ktorsion (int indSpecies=0)
	Allows the spring constant to be accessed. More...

void	set_rho (Mdouble new_, unsigned int indSpecies=0)
	Allows the density to be changed. More...

Mdouble	get_rho (int indSpecies=0)
	Allows the density to be accessed. More...

void	set_dispt (Mdouble new_, unsigned int indSpecies=0)
	Allows the tangential viscosity to be changed. More...

Mdouble	get_dispt (unsigned int indSpecies=0)
	Allows the tangential viscosity to be accessed. More...

void	set_disprolling (Mdouble new_, unsigned int indSpecies=0)
	Allows the tangential viscosity to be changed. More...

Mdouble	get_disprolling (unsigned int indSpecies=0)
	Allows the tangential viscosity to be accessed. More...

void	set_disptorsion (Mdouble new_, unsigned int indSpecies=0)
	Allows the tangential viscosity to be changed. More...

Mdouble	get_disptorsion (unsigned int indSpecies=0)
	Allows the tangential viscosity to be accessed. More...

void	set_disp (Mdouble new_, unsigned int indSpecies=0)
	Allows the normal dissipation to be changed. More...

Mdouble	get_disp (unsigned int indSpecies=0)
	Allows the normal dissipation to be accessed. More...

void	set_dissipation (Mdouble new_, unsigned int indSpecies=0)
	Allows the normal dissipation to be changed. More...

Mdouble	get_dissipation (unsigned int indSpecies=0)
	Allows the normal dissipation to be accessed. More...

void	set_mu (Mdouble new_, unsigned int indSpecies=0)
	Allows the Coulomb friction coefficient to be changed. More...

Mdouble	get_mu (unsigned int indSpecies=0)
	Allows the Coulomb friction coefficient to be accessed. More...

void	set_murolling (Mdouble new_, unsigned int indSpecies=0)
	Allows the Coulomb friction coefficient to be changed. More...

Mdouble	get_murolling (unsigned int indSpecies=0)
	Allows the Coulomb friction coefficient to be accessed. More...

void	set_mutorsion (Mdouble new_, unsigned int indSpecies=0)
	Allows the Coulomb friction coefficient to be changed. More...

Mdouble	get_mutorsion (unsigned int indSpecies=0)
	Allows the Coulomb friction coefficient to be accessed. More...

void	set_rotation (bool new_)

bool	get_rotation ()

void	set_dim_particle (int new_, unsigned int indSpecies=0)
	Allows the dimension of the particle (f.e. for mass) to be changed. More...

int	get_dim_particle (unsigned int indSpecies=0)
	Allows the dimension of the particle (f.e. for mass) to be accessed. More...

bool	get_save_data_data ()
	Returns the data counter. More...

bool	get_save_data_ene ()

bool	get_save_data_fstat ()

bool	get_save_data_stat ()

bool	get_do_stat_always ()

void	set_k_and_restitution_coefficient (Mdouble k_, Mdouble eps, Mdouble mass, unsigned int indSpecies=0)
	Sets k, disp such that it matches a given tc and eps for a collision of two copies of P. More...

void	set_collision_time_and_restitution_coefficient (Mdouble tc, Mdouble eps, Mdouble mass, unsigned int indSpecies=0)
	Sets k, disp such that it matches a given tc and eps for a collision of two copies of P. More...

void	set_collision_time_and_restitution_coefficient (Mdouble tc, Mdouble eps, Mdouble mass1, Mdouble mass2, unsigned int indSpecies=0)
	Set k, disp such that is matches a given tc and eps for a collision of two different masses. More...

void	set_collision_time_and_normal_and_tangential_restitution_coefficient (Mdouble tc, Mdouble eps, Mdouble beta, Mdouble mass1, Mdouble mass2, unsigned int indSpecies=0)
	See CSpecies::set_collision_time_and_normal_and_tangential_restitution_coefficient. More...

void	set_collision_time_and_normal_and_tangential_restitution_coefficient_nodispt (Mdouble tc, Mdouble eps, Mdouble beta, Mdouble mass1, Mdouble mass2, unsigned int indSpecies=0)
	See CSpecies::set_collision_time_and_normal_and_tangential_restitution_coefficient. More...

Mdouble	get_collision_time (Mdouble mass, unsigned int indSpecies=0)
	Calculates collision time for two copies of a particle of given disp, k, mass. More...

Mdouble	get_restitution_coefficient (Mdouble mass, unsigned int indSpecies=0)
	Calculates restitution coefficient for two copies of given disp, k, mass. More...

Mdouble	get_xmin ()
	Get xmin. More...

Mdouble	get_xmax ()
	Get xmax. More...

Mdouble	get_ymin ()
	Gets ymin. More...

Mdouble	get_ymax ()
	Gets ymax. More...

Mdouble	get_zmin ()
	Gets zmin. More...

Mdouble	get_zmax ()
	Gets zmax. More...

void	set_xmin (Mdouble new_xmin)
	Sets xmin and walls, assuming the standard definition of walls as in the default constructor. More...

void	set_ymin (Mdouble new_ymin)

void	set_zmin (Mdouble new_zmin)
	Sets ymin and walls, assuming the standard definition of walls as in the default constructor. More...

void	set_xmax (Mdouble new_xmax)
	Sets xmax and walls, assuming the standard definition of walls as in the default constructor. More...

void	set_ymax (Mdouble new_ymax)
	Sets ymax and walls, assuming the standard definition of walls as in the default constructor. More...

void	set_zmax (Mdouble new_zmax)
	Sets ymax and walls, assuming the standard definition of walls as in the default constructor. More...

void	set_dt (Mdouble new_dt)
	Allows the time step dt to be changed. More...

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

void	set_name (const char *name)
	Sets the name of the problem, used for the same data files. More...

void	set_xballs_colour_mode (int new_cmode)
	Set the xball output mode. More...

void	set_xballs_cmode (int new_cmode)

int	get_xballs_cmode ()

void	set_xballs_vector_scale (double new_vscale)
	Set the scale of vectors in xballs. More...

double	get_xballs_vscale ()

void	set_xballs_additional_arguments (std::string new_)
	Set the additional arguments for xballs. More...

std::string	get_xballs_additional_arguments ()

void	set_xballs_scale (Mdouble new_scale)
	Set the scale of the xballs problem. The default is fit to screen. More...

double	get_xballs_scale ()

void	set_gravity (Vec3D new_gravity)
	Allows the gravitational acceleration to be changed. More...

Vec3D	get_gravity ()
	Allows the gravitational acceleration to be accessed. More...

void	set_dim (int new_dim)
	Allows the dimension of the simulation to be changed. More...

int	get_dim ()
	Allows the dimension of the simulation to be accessed. More...

int	get_restart_version ()
	Gets restart_version. More...

void	set_restart_version (int new_)
	Sets restart_version. More...

bool	get_restarted ()
	Gets restarted. More...

Mdouble	get_max_radius ()
	Sets restarted. More...

void	set_restarted (bool new_)

bool	get_append ()
	Gets restarted. More...

void	set_append (bool new_)
	Sets restarted. More...

Mdouble	get_ene_ela ()
	Gets ene_ela. More...

void	set_ene_ela (Mdouble new_)
	Sets ene_ela. More...

void	add_ene_ela (Mdouble new_)
	Sets ene_ela. More...

void	Remove_Particle (int IP)
	This function removes partice IP from the vector of particles by moving the last particle in the vector to the position if IP Also it checks if the moved Particle has any tangentialsspring-information, which needs to be moved to a different particle, because tangential spring information always needs to be stored in the real particle with highest particle index. More...

Mdouble	get_Mass_from_Radius (Mdouble radius, int indSpecies=0)

Mdouble	get_maximum_velocity (BaseParticle &P)
	Calculates the maximum velocity allowed for a collision of two copies of P (for higher velocities particles could pass through each other) More...

virtual BaseParticle *	getSmallestParticle ()

virtual BaseParticle *	getLargestParticle ()

virtual void	removeParticle (int iP)

Mdouble	get_maximum_velocity ()

void	set_dt_by_mass (Mdouble mass)
	Sets dt to 1/50-th of the collision time for two particles of mass P. More...

void	set_dt (BaseParticle &P)
	Sets dt to 1/50-th of the collision time for two copies of P. More...

void	set_dt ()
	Sets dt to 1/50-th of the smallest possible collision time. More...

virtual void	setup_particles_initial_conditions ()
	This function allows the initial conditions of the particles to be set, by default locations is random. More...

virtual void	create_xballs_script ()
	This creates a scipt which can be used to load the xballs problem to display the data just generated. More...

virtual double	getInfo (BaseParticle &P)
	Allows the user to set what is written into the info column in the data file. By default is store the Species ID number. More...

virtual void	save_restart_data ()
	Stores all MD data. More...

int	load_restart_data ()
	Loads all MD data. More...

int	load_restart_data (std::string filename)

void	statistics_from_restart_data (const char *name)
	Loads all MD data and plots statistics for all timesteps in the .data file. More...

virtual void	write (std::ostream &os)
	Writes all MD data. More...

virtual void	read (std::istream &is)
	Reads all MD data. More...

virtual void	write_v1 (std::ostream &os)
	Writes all MD data. More...

virtual void	read_v1 (std::istream &is)
	Reads all MD data. More...

virtual void	read_v2 (std::istream &is)

bool	load_from_data_file (const char *filename, unsigned int format=0)
	This allows particle data to be reloaded from data files. More...

bool	load_par_ini_file (const char *filename)
	allows the user to read par.ini files (useful to read MDCLR files) More...

bool	read_next_from_data_file (unsigned int format=0)
	by default format do not pass an argument; only specify format if you have to read a special format (f.e. dim=2, but format=14 (3d format)) More...

int	read_dim_from_data_file ()

bool	find_next_data_file (Mdouble tmin, bool verbose=true)

virtual void	print (std::ostream &os, bool print_all=false)
	Outputs MD. More...

void	add_Species (CSpecies &S)

void	add_Species (void)

void	set_format (int new_)

int	get_format ()

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

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

Public Member Functions inherited from STD_save
	STD_save ()
	Default constructor: sets the counter to 0 (i.e. no number will be included). More...

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

void	constructor ()

void	inc_counter_in_file ()
	Increament the counter value stored in the file_counter by 1 and store the new value. More...

int	read_run_num_from_file ()
	Read rom the counter file the counter. More...

void	set_counter_from_file ()
	Sets the counter based on the current number stored in the counter file. More...

void	save_info_to_disk ()
	Saves the information generated by info to disk in a file. More...

void	set_counter (int new_counter)
	This set the counter, overriding the defaults. More...

int	get_counter ()
	This returns the current value of the counter. More...

bool	FileExists (std::string strFilename)
	Function to check if a file exists, is used to check if a run has already need done. More...

void	auto_number ()

std::vector< int >	get_numbers (int size_x, int size_y)
	This turns a counter into two indexs for doing parmater studies. The indexs run from 1:size_x and 1:size_y where as the study number starts at 0. More...

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

void	set_name (const char *name)
	Sets the name of the problem, used for the same data files. More...

std::string	get_name ()
	Allows the problem_name to be accessed. More...

std::fstream &	get_data_file ()
	Allows the problem_name to be accessed. More...

std::fstream &	get_stat_file ()
	Allows the problem_name to be accessed. More...

std::fstream &	get_fstat_file ()
	Allows the problem_name to be accessed. More...

std::fstream &	get_ene_file ()
	Allows the problem_name to be accessed. More...

void	set_fstat_filename (std::string filename)

void	set_data_filename (std::string filename)

void	set_stat_filename (std::string filename)

void	set_ene_filename (std::string filename)

void	set_fstat_filename ()

void	set_data_filename ()

void	set_stat_filename ()

void	set_ene_filename ()

std::string	get_fstat_filename ()

std::string	get_data_filename ()

std::string	get_stat_filename ()

std::string	get_ene_filename ()

void	set_step_size (unsigned int new_)

unsigned int	get_step_size ()

void	set_options_fstat (unsigned int new_)
	set and get for file options More...

unsigned int	get_options_fstat (void)

void	set_options_data (unsigned int new_)

unsigned int	get_options_data (void)

void	set_options_stat (unsigned int new_)

unsigned int	get_options_stat (void)

void	set_options_restart (unsigned int new_)

unsigned int	get_options_restart (void)

void	set_options_ene (unsigned int new_)

unsigned int	get_options_ene (void)

bool	open_file (std::fstream &file, std::string filename, unsigned int options, std::fstream::openmode mode)

bool	open_fstat_file (std::fstream::openmode mode=std::fstream::out)

bool	open_data_file (std::fstream::openmode mode=std::fstream::out)

bool	open_stat_file (std::fstream::openmode mode=std::fstream::out)

bool	open_ene_file (std::fstream::openmode mode=std::fstream::out)

bool	open_counted_file (std::fstream &file, std::string filenameNoCount, std::fstream::openmode mode)
	opens file needed if data is written in multiple files More...

bool	increase_counter_fstat (std::fstream::openmode mode)

bool	increase_counter_data (std::fstream::openmode mode)

bool	increase_counter_stat (std::fstream::openmode mode)

bool	increase_counter_ene (std::fstream::openmode mode)

void	set_file_counter (int new_)

int	get_file_counter ()

Protected Member Functions
bool	satisfiesInclusionCriteria (BaseParticle *P)

bool	loadVelocityProfile (const char *filename)

Vec3D	getVelocityProfile (Vec3D Position)

bool	read_next_from_p3p_file ()

void	auto_setdim ()

Protected Member Functions inherited from MD
virtual void	compute_all_forces ()
	This does the force computation. More...

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

CTangentialSpring *	getTangentialSpring (BaseParticle PI, BaseParticle PJ, BaseParticle *PJreal)

CTangentialSpring *	getTangentialSpringWall (BaseParticle *pI, int w)

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

void	compute_plastic_internal_forces (BaseParticle P1, BaseParticle P2)
	Computes plastic forces between particles. More...

virtual void	compute_external_forces (BaseParticle *PI)
	This is were the computation of external forces takes place (e.g. gravity) More...

virtual void	compute_walls (BaseParticle *PI)
	This is were the walls are. More...

Mdouble	computeShortRangeForceWithWall (BaseParticle pI, int wall, CSpecies pSpecies, Mdouble dist)

Mdouble	computeShortRangeForceWithParticle (BaseParticle PI, BaseParticle PJ, BaseParticle PJreal, CSpecies pSpecies, Mdouble dist)

virtual void	actions_before_time_loop ()
	This is actions before the start of the main time loop. More...

virtual void	HGRID_actions_before_time_loop ()
	This is actions before the start of the main time loop. More...

virtual void	HGRID_actions_before_time_step ()
	This is action before the time step is started. More...

virtual void	HGRID_InsertParticleToHgrid (BaseParticle *obj UNUSED)
	This is action before the time step is started. More...

virtual void	HGRID_UpdateParticleInHgrid (BaseParticle *obj UNUSED)

virtual void	HGRID_RemoveParticleFromHgrid (BaseParticle *obj UNUSED)

virtual bool	get_HGRID_UpdateEachTimeStep ()

virtual void	actions_before_time_step ()
	This is action before the time step is started. More...

virtual void	actions_after_solve ()
	This is actions at the end of the code, but before the files are closed. More...

virtual void	actions_after_time_step ()
	This is action after the time step is finished. More...

virtual void	output_xballs_data ()
	Output xball data for Particle i. More...

virtual void	output_xballs_data_particle (int i)
	This function outputs the location and velocity of the particle in a format the xballs progream can read. More...

virtual void	start_ene ()
	Functions for ene file. More...

virtual void	fstat_header ()

virtual void	output_ene ()
	This function outputs statistical data - The default is to compute the rotational kinetic engergy, linear kinetic energy, and the centre of mass. More...

virtual void	set_initial_pressures_for_pressure_controlled_walls ()

virtual void	do_integration_before_force_computation (BaseParticle *pI)
	This is were the integration is done. More...

virtual void	checkInteractionWithBoundaries ()

virtual void	HGRID_update_move (BaseParticle *, Mdouble)

virtual void	HGRID_actions_before_integration ()

virtual void	HGRID_actions_after_integration ()

virtual void	do_integration_after_force_computation (BaseParticle *pI)
	This is were the integration is done. More...

virtual void	InitBroadPhase ()
	Initialisation of Broad Phase Information (Default no Broad Phase so empty) More...

virtual void	broad_phase (BaseParticle *i)
	Broad phase of contact detection goes here. Default check all contacts. More...

void	set_FixedParticles (unsigned int n)

void	initialize_tangential_springs ()

void	compute_particle_masses ()
	Computes the mass of each particle. More...

virtual void	cout_time ()
	std::couts time More...

virtual bool	continue_solve ()

void	reset_DeltaMax ()
	sets the history parameter DeltaMax of all particles to zero More...

void	reset_TangentialSprings ()
	sets the history parameter TangentialSprings of all particles to zero 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 points of evaluation are placed in an grid on the domain [xminStat,xmaxStat]x[yminStat,ymaxStat]x[zminStat,zmaxStat]. 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
	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 More...

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

int	format

Mdouble	mirrorAtDomainBoundary

bool	isMDCLR

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

bool	doublePoints

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

Protected Attributes inherited from MD
std::vector< CSpecies >	Species
	These are the particle parameters like dissipation etc. More...

Protected Attributes inherited from STD_save
std::stringstream	problem_name
	Stores the problem_name. More...

std::stringstream	data_filename
	These store the save file names, by default they are derived from problem_name. More...

std::stringstream	stat_filename

std::stringstream	fstat_filename

std::stringstream	ene_filename

std::fstream	data_file
	Stream used for data files. More...

std::fstream	stat_file

std::fstream	fstat_file

std::fstream	ene_file

unsigned int	options_fstat
	Indicators if files are created or not 0: file will not be created 1: file will be written in one file 2: file will be written in multiple files. More...

unsigned int	options_data

unsigned int	options_stat

unsigned int	options_ene

unsigned int	options_restart

unsigned int	file_counter
	Counter needed if file will be written in multiple files. More...

unsigned int	step_size

Additional Inherited Members
Public Attributes inherited from MD
RNG	random

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 30 of file StatisticsPoint.h.

Constructor & Destructor Documentation

template<StatType T>

StatisticsVector< T >::StatisticsVector ( )

inline

Basic constructor only calls constructor()

Definition at line 56 of file StatisticsVector.h.

References StatisticsVector< T >::constructor().

56 {constructor();}

StatisticsVector::constructor

void constructor()

Sets up all basic things.

Definition: StatisticsVector.hcc:74

template<StatType T>

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

inline

Definition at line 57 of file StatisticsVector.h.

References StatisticsVector< T >::constructor().

57 {constructor(name);}

StatisticsVector::constructor

void constructor()

Sets up all basic things.

Definition: StatisticsVector.hcc:74

template<StatType T>

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

Copy constructor.

Definition at line 139 of file StatisticsVector.hcc.

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

template<StatType T>

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

Advanced constructor that accepts arguments from the command line.

Definition at line 224 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 948 of file StatisticsVector.hcc.

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

template<StatType T>

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

Output average of statistical variables.

Output names of statistical variables.

Definition at line 649 of file StatisticsVector.hcc.

References StatisticsPoint< T >::set_zero().

                                                                                  {
     StatisticsPoint<T> avg;
     avg.set_zero();
     for (unsigned int i=0; i<P.size(); i++) avg += P[i];
     avg /= P.size();
     return avg;
 }

template<StatType T>

bool StatisticsVector< T >::check_current_time_for_statistics ( )

inline

Definition at line 91 of file StatisticsVector.h.

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

91 {return (get_t()>get_tminStat()&&get_t()<=get_tmaxStat()+get_dt());};

StatisticsVector::get_tminStat

Mdouble get_tminStat()

Definition: StatisticsVector.h:88

MD::get_dt

Mdouble get_dt()

Allows the time step dt to be accessed.

Definition: MD.h:339

MD::get_t

Mdouble get_t()

Access function for the time.

Definition: MD.h:108

StatisticsVector::get_tmaxStat

Mdouble get_tmaxStat()

Definition: StatisticsVector.h:89

template<StatType T>

void StatisticsVector< T >::constructor ( )

Sets up all basic things.

Definition at line 74 of file StatisticsVector.hcc.

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

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

                                       { 
     //stores the template parameter Stattype in a variable
     set_statType();
     
     // sets connection between StatisticsPoint and StatisticsVector
     StatisticsPoint<T>::set_gb(this);
     // set nx, ny,nz
     nx = ny = nz = 1;
     nxMirrored=nyMirrored=nzMirrored=0;
 
     // set default CG variables 
     CG_type = Gaussian;
     // in get_statistics_..., w2 is set to w_over_rmax*rmax if w2==0
     w_over_rmax = 1;
     
     // bounded domain
     //~ boundedDomain = false;
     
     // set default values; variables will only be used if the default 
     // is overwritten; thus the default is a nonsense value 
     w2 = 0.0;
     set_tintStat(0);
     set_tminStat(-1e20);
     set_tmaxStat(NAN);
     rmin = 0;
     hmax = 1e20;
     rmax = 1e20;
     indSpecies = -1; //all species
 
     
     //set default options - Not used it physially don't make sense i.e if ?max<?min
     xminStat = NAN;
     xmaxStat = NAN;
     yminStat = NAN;
     ymaxStat = NAN;
     zminStat = NAN;
     zmaxStat = NAN; 
     
 
     // set default options
     ignoreFixedParticles = true; //this should be true if you want statistics only of the flow
     StressTypeForFixedParticles = 1;
     verbosity = 1;
     walls = true;
     periodicWalls = true;
     format = 0;
     mirrorAtDomainBoundary = false;
     set_superexact(false);
     VelocityProfile.resize(0); //equivalent to do not use velocityprofile
 
     // time averaging
     set_doTimeAverage(true);
     nTimeAverage = 0;
     nTimeAverageReset = -1;
     //calculate variance
     set_doVariance(false);
     //calculate gradient
     set_doGradient(false);
     doublePoints=false;
 
     // additional stuff
     set_options_stat(1);
 }

template<StatType T>

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

Definition at line 202 of file StatisticsVector.hcc.

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

template<StatType T>

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

get force statistics

todo{What are P1 and P2 over here?}

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

Definition at line 1719 of file StatisticsVector.hcc.

References Cross(), PeriodicBoundary::distance(), PeriodicBoundary::shift_positions(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

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

template<StatType T>

void StatisticsVector< T >::evaluate_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}

Definition at line 1851 of file StatisticsVector.hcc.

References Cross(), PeriodicBoundary::distance(), Dyadic(), MD::get_t(), PeriodicBoundary::shift_position(), sqr, SymmetrizedDyadic(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

                                                                                                  {
     bool doDisplacement = false;
     
     if (periodicWalls)
         for (unsigned int k=wp; k<getBoundaryHandler().getNumberOfObjects(); k++)
         {
             PeriodicBoundary* b0=dynamic_cast<PeriodicBoundary*>(getBoundaryHandler().getObject(k));
             if(b0)
             {
                 Mdouble dist = b0->distance(P1);
                 if (sqr(dist-(*P)->get_Radius())<9.0*get_w2())
                 {
                     b0->shift_position(*P);
                     evaluate_particle_statistics(P, k+1);
                     b0->shift_position(*P);
                 }               
             }
         }
 
     //evaluate fields that only depend on CParticle parameters
     Mdouble phi; //Course graining function
     Vec3D V = Vec3D(0,0,0);
     if (VelocityProfile.size()) V = getVelocityProfile ((*P)->get_Position());
 
     for (unsigned int i=0; i<Points.size(); i++) {
         phi = Points[i].CG_function((*P)->get_Position());
         if (phi) {
             if (!std::isfinite(phi)) {
                 std::cout 
                 << "t =" << MD::get_t()
                 << "error: phi =" << phi 
                 << " is infinite for P=" << (*P)->get_Position()
                 << ", Point=" << Points[i].get_Position()
                 << std::endl;
                 phi=0;
             }
             Points[i].Nu += (*P)->get_Volume(get_Species()) * phi;
             Points[i].Density += (*P)->get_Mass() * phi;
             if (VelocityProfile.size()) {
                 Points[i].Momentum += ((*P)->get_Velocity()-V) * ((*P)->get_Mass() * phi);
                 Points[i].MomentumFlux += SymmetrizedDyadic((*P)->get_Velocity()-V, (*P)->get_Velocity()-V) * ((*P)->get_Mass() * phi);
             } else {
                 Points[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * phi);
                 Points[i].MomentumFlux += SymmetrizedDyadic((*P)->get_Velocity(), (*P)->get_Velocity()) * ((*P)->get_Mass() * phi);
             }
             Points[i].DisplacementMomentum += (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()) * ((*P)->get_Mass() * phi);
             Points[i].DisplacementMomentumFlux += SymmetrizedDyadic((*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()), (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount())) * ((*P)->get_Mass() * phi);
             Points[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * phi);
             
             Vec3D LocalAngularMomentum = Cross(Points[0].clearAveragedDirections((*P)->get_Position()-Points[i].get_Position()), (*P)->get_Velocity()) * (*P)->get_Mass() 
                 + (*P)->get_AngularVelocity() * (*P)->get_Inertia();
             Points[i].LocalAngularMomentum += phi * LocalAngularMomentum;
             Points[i].LocalAngularMomentumFlux += Dyadic(LocalAngularMomentum, (*P)->get_Velocity() * phi);
 
             //Note: Displacement is only computed directly if gradients are cmputed
             if (get_doGradient()) {
                 
                 Vec3D d_phi = Points[i].CG_gradient((*P)->get_Position(), phi);
                 
                 if (doDisplacement) {
                     //begin: Linear displacement, \f$1/(2 \rho^2) \sum_{pq} m_p m_q \phi_q *(v_{pqa} \partial_b \phi_p + v_{pqb} \partial_a \phi_p) \f$ 
                     for (std::vector<BaseParticle*>::iterator Q = getParticleHandler().begin(); Q!=getParticleHandler().end(); ++Q) {
                         double phiQ = Points[i].CG_function((*Q)->get_Position()); //Course graining function
                         Points[i].Displacement += SymmetrizedDyadic((*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()) - (*Q)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()), d_phi) 
                          * ((*P)->get_Mass() * (*Q)->get_Mass() * phiQ);
                     }
                     //end:Displacement
                 }
                 
                 dx[i].Nu += (*P)->get_Volume(get_Species()) * d_phi.X;
                 dx[i].Density += (*P)->get_Mass() * d_phi.X;
                 dx[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * d_phi.X);
                 dx[i].DisplacementMomentum += (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()) * ((*P)->get_Mass() * d_phi.X);
                 dx[i].MomentumFlux += SymmetrizedDyadic((*P)->get_Velocity(), (*P)->get_Velocity()) * ((*P)->get_Mass() * d_phi.X);
                 dx[i].DisplacementMomentumFlux += SymmetrizedDyadic((*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()), (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount())) * ((*P)->get_Mass() * d_phi.X);
                 dx[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * d_phi.X);
 
                 dy[i].Nu += (*P)->get_Volume(get_Species()) * d_phi.Y;
                 dy[i].Density += (*P)->get_Mass() * d_phi.Y;
                 dy[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * d_phi.Y);
                 dy[i].DisplacementMomentum += (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()) * ((*P)->get_Mass() * d_phi.Y);
                 dy[i].MomentumFlux += SymmetrizedDyadic((*P)->get_Velocity(), (*P)->get_Velocity()) * ((*P)->get_Mass() * d_phi.Y);
                 dy[i].DisplacementMomentumFlux += SymmetrizedDyadic((*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()), (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount())) * ((*P)->get_Mass() * d_phi.Y);
                 dy[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * d_phi.Y);
 
                 dz[i].Nu += (*P)->get_Volume(get_Species()) * d_phi.Z;
                 dz[i].Density += (*P)->get_Mass() * d_phi.Z;
                 dz[i].Momentum += (*P)->get_Velocity() * ((*P)->get_Mass() * d_phi.Z);
                 dz[i].DisplacementMomentum += (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()) * ((*P)->get_Mass() * d_phi.Z);
                 dz[i].MomentumFlux += SymmetrizedDyadic((*P)->get_Velocity(), (*P)->get_Velocity()) * ((*P)->get_Mass() * d_phi.Z);
                 dz[i].DisplacementMomentumFlux += SymmetrizedDyadic((*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount()), (*P)->get_Displacement2(get_xmin(),get_xmax(),get_ymin(),get_ymax(),get_zmin(),get_zmax(),get_dt()*get_savecount())) * ((*P)->get_Mass() * d_phi.Z);
                 dz[i].EnergyFlux += (*P)->get_Velocity() * ((*P)->get_Mass()*(*P)->get_Velocity().GetLength2()/2 * d_phi.Z);
             } //end if get_doGradient
         } //end if phi
     } // end forall Points
 }

template<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 1787 of file StatisticsVector.hcc.

References PeriodicBoundary::distance(), PeriodicBoundary::shift_positions(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

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

template<StatType T>

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

inline

Definition at line 252 of file StatisticsVector.h.

References StatisticsVector< T >::CGPolynomial.

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

template<StatType T>

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

inline

Definition at line 244 of file StatisticsVector.h.

References StatisticsVector< T >::CGPolynomial.

                                           {
         return CGPolynomial.evaluate(r);
     }

template<StatType T>

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

inline

Definition at line 248 of file StatisticsVector.h.

References StatisticsVector< T >::CGPolynomial.

                                                   {
         return CGPolynomial.evaluateGradient(r);
     }

template<StatType T>

void StatisticsVector< T >::finish_statistics ( )

virtual

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

Reimplemented from MD.

Definition at line 729 of file StatisticsVector.hcc.

 {
     if (get_doTimeAverage()) write_time_average_statistics();
     stat_file.close();
 }

template<StatType T>

void StatisticsVector< T >::gather_force_statistics_from_fstat_and_data ( )

get force statistics from particle collisions

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

Definition at line 1662 of file StatisticsVector.hcc.

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

template<StatType T>

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

get force statistics from particle collisions

Todo:: This has to be tested

Definition at line 1172 of file StatisticsVector.hcc.

References BaseParticle::get_Position(), and BaseParticle::get_Radius().

 {
     std::string dummyStr;
     Mdouble N, dummy;
     
     //Ignoring first line
     getline(p3c_file, dummyStr); 
     std::cout << dummyStr << std::endl; 
     //Reading second line
     p3c_file >> dummy; set_t(dummy);
     p3c_file >> N;
     std::cout << "t= " << get_t() << ": loading " << N << " contacts ..." << std::endl;
     getline(p3c_file,dummyStr);
     //Ignoring third line
     getline(p3c_file,dummyStr);
     std::cout << dummyStr << std::endl;
     //Checking for end of file
     //if (p3c_file.eof()||p3c_file.peek()==-1) return false;
     
     //set up index list
     int Nmax = 0;
     for (std::vector<BaseParticle*>::iterator it = getParticleHandler().begin(); it!=getParticleHandler().end(); ++it) 
         Nmax = std::max(Nmax,(*it)->get_Index());
     static std::vector<int> index;
     index.resize(Nmax);
     int i=0;
     for (std::vector<BaseParticle*>::iterator it = getParticleHandler().begin(); it!=getParticleHandler().end(); ++it) {
         index[(*it)->get_Index()]=i;
         i++;
     }
     
     // ==> Particledata_4.2to4.25s.csv.p4c <==
     // TIMESTEP CONTACTS
     // 4.2001  31367
     //  P1 P2 CX CY CZ FX FY FZ
     //  9711 9424 0.0549396 0.009705380.149582 4.37283e-05 2.02347e-05 -0.000250662
 
     // ==> /Users/weinhartt/CarlosLabra/ConfinedCompression_P3P_Spheres.csv.p3c <==
     // TIMESTEP CONTACTS
     // 6.500010000000e-01  21886
     //  P1 P2 FX FY FZ
     //  1126 218 -1.899612000000e-02 6.716474700000e-02 -1.527340000000e-02
 
     int index1_p3, index2_p3;
     int index1, index2;
     Vec3D Contact, Force;
     Mdouble delta, fdotn, fdott;
     Vec3D P1_P2_normal, P1_P2_tangential;
     Mdouble dist;
     BaseParticle* P1;
     BaseParticle* P2;
 
     Vec3D contact;
     int counter = 0;
     //go through N lines in the p4c file
     for (int i=0; i<N; i++) {
         counter++;
         if (version==3) {
             p3c_file 
                 >> index1_p3
                 >> index2_p3
                 >> Force;
         } else {
             p3c_file 
                 >> index1_p3
                 >> index2_p3
                 >> Contact
                 >> Force;
         }
         index1=index[index1_p3];
         index2=index[index2_p3];
         P1 = getParticleHandler().getObject(index1);
         P2 = getParticleHandler().getObject(index2);
         P1_P2_normal = P1->get_Position() - P2->get_Position();
         dist = GetLength(P1_P2_normal);
         P1_P2_normal /= dist;
         delta = P1->get_Radius() + P2->get_Radius() - dist;
         if (version==3) {
             Contact = P2->get_Position() + P1_P2_normal * (P2->get_Radius()-delta/2);
         }
         fdotn = -Dot(Force,P1_P2_normal)/GetLength2(P1_P2_normal);
         P1_P2_tangential = Force + fdotn * P1_P2_normal; //get tangential force
         fdott = GetLength(P1_P2_tangential);
         if (fdott==0) {
             P1_P2_tangential = Vec3D(0,0,0);
         } else  {   
             P1_P2_tangential /= fdott;
         }
         //Finished reading file
         gather_statistics_collision(index1,index2,Contact,delta,0,fdotn,fdott, P1_P2_normal,-P1_P2_tangential);
         gather_statistics_collision(index2,index1,Contact,delta,0,fdotn,fdott, -P1_P2_normal, P1_P2_tangential);
     }
     
     
     if (verbosity>1) std::cout << "#forces="<< counter << std::endl;
     //Ignoring last line
     getline(p3c_file,dummyStr);
     std::cout << dummyStr << std::endl;
     
     gather_force_statistics_from_p3w(version,index);
 }

template<StatType T>

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

get force statistics from particle collisions

Todo:: This has to be tested

Definition at line 1276 of file StatisticsVector.hcc.

References BaseParticle::get_Position(), and BaseParticle::get_Radius().

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

template<StatType T>

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

virtual

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

Todo:: : this is to make the stresses infinitely (i.e. 300 units) long; has to be improved 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 MD.

Definition at line 1518 of file StatisticsVector.hcc.

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

template<StatType T>

CG StatisticsVector< T >::get_CG_type ( )

inline

Definition at line 95 of file StatisticsVector.h.

References StatisticsVector< T >::CG_type.

95 {return CG_type;};

StatisticsVector::CG_type

CG CG_type

coarse graining type (Gaussian, Heaviside, Polynomial)

Definition: StatisticsVector.h:309

template<StatType T>

Mdouble StatisticsVector< T >::get_cutoff ( )

inline

Definition at line 108 of file StatisticsVector.h.

References StatisticsVector< T >::cutoff.

108 {return cutoff;};

StatisticsVector::cutoff

Mdouble cutoff

The distance from the origin at which the cg function vanishes; cutoff=w for HeavisideSphere or Polyn...

Definition: StatisticsVector.h:317

template<StatType T>

Mdouble StatisticsVector< T >::get_cutoff2 ( )

inline

Definition at line 109 of file StatisticsVector.h.

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

109 {return sqr(cutoff);};

StatisticsVector::cutoff

Mdouble cutoff

The distance from the origin at which the cg function vanishes; cutoff=w for HeavisideSphere or Polyn...

Definition: StatisticsVector.h:317

sqr

#define sqr(a)

Definition: ExtendedMath.h:36

template<StatType T>

bool StatisticsVector< T >::get_doGradient ( )

inline

Definition at line 144 of file StatisticsVector.h.

References StatisticsVector< T >::doGradient.

144 {return doGradient;};

StatisticsVector::doGradient

bool doGradient

Determines if gradient is outputted.

Definition: StatisticsVector.h:303

template<StatType T>

bool StatisticsVector< T >::get_doTimeAverage ( )

inline

Definition at line 129 of file StatisticsVector.h.

References StatisticsVector< T >::doTimeAverage.

129 {return doTimeAverage;};

StatisticsVector::doTimeAverage

bool doTimeAverage

Determines if output is averaged over time.

Definition: StatisticsVector.h:297

template<StatType T>

bool StatisticsVector< T >::get_doublePoints ( )

inline

Definition at line 234 of file StatisticsVector.h.

References StatisticsVector< T >::doublePoints.

234 {return doublePoints;}

StatisticsVector::doublePoints

bool doublePoints

Definition: StatisticsVector.h:363

template<StatType T>

bool StatisticsVector< T >::get_doVariance ( )

inline

Definition at line 141 of file StatisticsVector.h.

References StatisticsVector< T >::doVariance.

141 {return doVariance;};

StatisticsVector::doVariance

bool doVariance

Determines if variance is outputted.

Definition: StatisticsVector.h:301

template<StatType T>

bool StatisticsVector< T >::get_ignoreFixedParticles ( )

inline

Definition at line 150 of file StatisticsVector.h.

References StatisticsVector< T >::ignoreFixedParticles.

150 {return ignoreFixedParticles;};

StatisticsVector::ignoreFixedParticles

bool ignoreFixedParticles

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

Definition: StatisticsVector.h:344

template<StatType T>

Mdouble StatisticsVector< T >::get_mirrorAtDomainBoundary ( )

inline

Definition at line 138 of file StatisticsVector.h.

References StatisticsVector< T >::mirrorAtDomainBoundary.

138 {return mirrorAtDomainBoundary;};

StatisticsVector::mirrorAtDomainBoundary

Mdouble mirrorAtDomainBoundary

Definition: StatisticsVector.h:356

template<StatType T>

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

inline

Definition at line 98 of file StatisticsVector.h.

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

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

StatisticsVector::nx

int nx

Grid size nx,ny,nz (by default the points of evaluation are placed in an grid on the domain [xminStat...

Definition: StatisticsVector.h:264

StatisticsVector::ny

int ny

see nx

Definition: StatisticsVector.h:266

StatisticsVector::nz

int nz

see nx

Definition: StatisticsVector.h:268

template<StatType T>

int StatisticsVector< T >::get_nTimeAverage ( )

inline

Definition at line 212 of file StatisticsVector.h.

References StatisticsVector< T >::nTimeAverage.

212 {return nTimeAverage;};

StatisticsVector::nTimeAverage

int nTimeAverage

A counter needed to average over time.

Definition: StatisticsVector.h:305

template<StatType T>

int StatisticsVector< T >::get_nx ( )

inline

Definition at line 82 of file StatisticsVector.h.

References StatisticsVector< T >::nx.

82 {return nx;};

StatisticsVector::nx

int nx

Grid size nx,ny,nz (by default the points of evaluation are placed in an grid on the domain [xminStat...

Definition: StatisticsVector.h:264

template<StatType T>

int StatisticsVector< T >::get_ny ( )

inline

Definition at line 83 of file StatisticsVector.h.

References StatisticsVector< T >::ny.

83 {return ny;};

StatisticsVector::ny

int ny

see nx

Definition: StatisticsVector.h:266

template<StatType T>

int StatisticsVector< T >::get_nz ( )

inline

Definition at line 84 of file StatisticsVector.h.

References StatisticsVector< T >::nz.

84 {return nz;};

StatisticsVector::nz

int nz

see nx

Definition: StatisticsVector.h:268

template<StatType T>

bool StatisticsVector< T >::get_periodicWalls ( )

inline

Definition at line 160 of file StatisticsVector.h.

References StatisticsVector< T >::periodicWalls.

160 {return periodicWalls;};

StatisticsVector::periodicWalls

bool periodicWalls

Turns off periodic walls before evaluation (needed for averaging, because we do not yet check if part...

Definition: StatisticsVector.h:342

template<StatType T>

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

inline

Definition at line 197 of file StatisticsVector.h.

References StatisticsVector< T >::Points.

197 {return Points;};

StatisticsVector::Points

std::vector< StatisticsPoint< T > > Points

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

Definition: StatisticsVector.h:277

template<StatType T>

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

inline

Definition at line 228 of file StatisticsVector.h.

References StatisticsVector< T >::CGPolynomial.

                                    {    
         return CGPolynomial.get_name();
     }

template<StatType T>

int StatisticsVector< T >::get_StressTypeForFixedParticles ( )

inline

Definition at line 132 of file StatisticsVector.h.

References StatisticsVector< T >::StressTypeForFixedParticles.

132 {return StressTypeForFixedParticles;};

StatisticsVector::StressTypeForFixedParticles

int StressTypeForFixedParticles

0 no Stress from fixed particles 1 Stress from fixed particles distributed between Contact and flowin...

Definition: StatisticsVector.h:349

template<StatType T>

bool StatisticsVector< T >::get_superexact ( )

inline

Definition at line 147 of file StatisticsVector.h.

References StatisticsVector< T >::superexact.

147 {return superexact;};

StatisticsVector::superexact

bool superexact

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

Definition: StatisticsVector.h:360

template<StatType T>

bool StatisticsVector< T >::get_TimeAverageReset ( )

inline

Definition at line 238 of file StatisticsVector.h.

References StatisticsVector< T >::nTimeAverageReset.

238 {return nTimeAverageReset;}

StatisticsVector::nTimeAverageReset

int nTimeAverageReset

Determines after how many timesteps the time average is reset.

Definition: StatisticsVector.h:299

template<StatType T>

Mdouble StatisticsVector< T >::get_tintStat ( )

inline

Definition at line 90 of file StatisticsVector.h.

References StatisticsVector< T >::tintStat.

90 {return tintStat;};

StatisticsVector::tintStat

Mdouble tintStat

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

Definition: StatisticsVector.h:328

template<StatType T>

Mdouble StatisticsVector< T >::get_tmaxStat ( )

inline

Definition at line 89 of file StatisticsVector.h.

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

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

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

StatisticsVector::tmaxStat

Mdouble tmaxStat

Statistical output will only be created if t

Definition: StatisticsVector.h:326

MD::get_tmax

Mdouble get_tmax()

Allows the upper time limit to be accessed.

Definition: MD.h:144

template<StatType T>

Mdouble StatisticsVector< T >::get_tminStat ( )

inline

Definition at line 88 of file StatisticsVector.h.

References StatisticsVector< T >::tminStat.

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

88 {return tminStat;};

StatisticsVector::tminStat

Mdouble tminStat

Statistical output will only be created if t>tminStat.

Definition: StatisticsVector.h:324

template<StatType T>

int StatisticsVector< T >::get_verbosity ( )

inline

Definition at line 154 of file StatisticsVector.h.

References StatisticsVector< T >::verbosity.

154 {return verbosity;};

StatisticsVector::verbosity

int verbosity

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

Definition: StatisticsVector.h:354

template<StatType T>

Mdouble StatisticsVector< T >::get_w ( )

inline

Definition at line 106 of file StatisticsVector.h.

References StatisticsVector< T >::w2.

106 {return sqrt(w2);};

StatisticsVector::w2

Mdouble w2

coarse graining width squared; for HeavisideSphere and Gaussian

Definition: StatisticsVector.h:315

template<StatType T>

Mdouble StatisticsVector< T >::get_w2 ( )

inline

Definition at line 107 of file StatisticsVector.h.

References StatisticsVector< T >::w2.

107 {return w2;};

StatisticsVector::w2

Mdouble w2

coarse graining width squared; for HeavisideSphere and Gaussian

Definition: StatisticsVector.h:315

template<StatType T>

Mdouble StatisticsVector< T >::get_w_over_rmax ( )

inline

Definition at line 163 of file StatisticsVector.h.

References StatisticsVector< T >::w_over_rmax.

163 {return w_over_rmax;};

StatisticsVector::w_over_rmax

Mdouble w_over_rmax

if w is not set manually then w will be set by multiplying this value by the largest particle radius ...

Definition: StatisticsVector.h:319

template<StatType T>

bool StatisticsVector< T >::get_walls ( )

inline

Definition at line 157 of file StatisticsVector.h.

References StatisticsVector< T >::walls.

157 {return walls;};

StatisticsVector::walls

bool walls

Turns off walls before evaluation.

Definition: StatisticsVector.h:339

template<StatType T>

Mdouble StatisticsVector< T >::get_xmaxStat ( )

inline

Definition at line 203 of file StatisticsVector.h.

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

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

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

StatisticsVector::xmaxStat

Mdouble xmaxStat

Definition: StatisticsVector.h:271

MD::get_xmax

Mdouble get_xmax()

Get xmax.

Definition: MD.h:307

template<StatType T>

Mdouble StatisticsVector< T >::get_xminStat ( )

inline

Functions to acces and change the domain of statistics.

Definition at line 200 of file StatisticsVector.h.

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

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

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

StatisticsVector::xminStat

Mdouble xminStat

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

Definition: StatisticsVector.h:271

MD::get_xmin

Mdouble get_xmin()

Get xmin.

Definition: MD.h:305

template<StatType T>

Mdouble StatisticsVector< T >::get_ymaxStat ( )

inline

Definition at line 204 of file StatisticsVector.h.

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

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

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

MD::get_ymax

Mdouble get_ymax()

Gets ymax.

Definition: MD.h:311

StatisticsVector::ymaxStat

Mdouble ymaxStat

Definition: StatisticsVector.h:271

template<StatType T>

Mdouble StatisticsVector< T >::get_yminStat ( )

inline

Definition at line 201 of file StatisticsVector.h.

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

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

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

StatisticsVector::yminStat

Mdouble yminStat

Definition: StatisticsVector.h:271

MD::get_ymin

Mdouble get_ymin()

Gets ymin.

Definition: MD.h:309

template<StatType T>

Mdouble StatisticsVector< T >::get_zmaxStat ( )

inline

Definition at line 205 of file StatisticsVector.h.

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

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

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

StatisticsVector::zmaxStat

Mdouble zmaxStat

Definition: StatisticsVector.h:271

MD::get_zmax

Mdouble get_zmax()

Gets zmax.

Definition: MD.h:315

template<StatType T>

Mdouble StatisticsVector< T >::get_zminStat ( )

inline

Definition at line 202 of file StatisticsVector.h.

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

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

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

MD::get_zmin

Mdouble get_zmin()

Gets zmin.

Definition: MD.h:313

StatisticsVector::zminStat

Mdouble zminStat

Definition: StatisticsVector.h:271

template<StatType T>

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

protected

Definition at line 307 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 ((double)(Position.X-VelocityProfile_Min.X)/VelocityProfile_D.X, &i);  
     m = modf ((double)(Position.Y-VelocityProfile_Min.Y)/VelocityProfile_D.Y, &j);  
     n = modf ((double)(Position.Z-VelocityProfile_Min.Z)/VelocityProfile_D.Z, &k);
     if (nx<1) i=0;
     if (ny<1) j=0;
     if (ny<1) k=0;
     int ix = (i*ny+j)*nz+k;
     Vec3D V;
     if (nx<=1) V.X = VelocityProfile[ix].X;
     else V.X=VelocityProfile[ix].X*(1-l)+VelocityProfile[ix+ny*nz].X*l;
     if (ny<=1) V.Y = VelocityProfile[ix].Y;
     else V.Y=VelocityProfile[ix].Y*(1-m)+VelocityProfile[ix+nz].Y*m;
     if (nz<=1) V.Z = VelocityProfile[ix].Z;
     else V.Z=VelocityProfile[ix].Z*(1-n)+VelocityProfile[ix+1].Z*n;
     return V;
 }

template<StatType T>

void StatisticsVector< T >::initialize_statistics ( )

virtual

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

This creates the file statistics will be saved to

Reimplemented from MD.

Definition at line 658 of file StatisticsVector.hcc.

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

 {
     reset_statistics();
 
     StatisticsVector<T>::set_w2(StatisticsVector<T>::get_w2());
         
     StatisticsVector<T>::set_Positions();
     
     //set tmin and tmax if tint is set
     if (get_tintStat()) {
         set_tminStat(get_t());
         set_tmaxStat(get_tminStat() + get_tintStat());
         std::cout << "tint set, time goes from: " << get_tminStat() << " to: " << get_tmaxStat() << std::endl;
     }
 
     if (!strcmp(get_stat_filename().c_str(),"")) set_stat_filename(); //set ouput filename to default unless otherwise specified
     if(!open_stat_file(std::fstream::out)) {std::cerr << "Problem opening stat file aborting" <<std::endl; exit(-1);};
         
     get_stat_file() << Points.begin()->write_variable_names() << std::endl;
     get_stat_file() << print_CG() << std::endl;
     
     //std::couts variables
     if (verbosity>1) std::cout << std::endl;
     if (verbosity>0) std::cout << std::endl << print() << std::endl;
 
 }

template<StatType T>

void StatisticsVector< T >::jump_fstat ( )

Definition at line 1493 of file StatisticsVector.hcc.

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

template<StatType T>

void StatisticsVector< T >::jump_p3c ( )

get force statistics from particle collisions

Definition at line 1136 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; set_t(dummy);
     p3c_file >> N;
     std::cout << "ignoring t= " << get_t() << ": loading " << N << " contacts ..." << std::endl;
     getline(p3c_file,dummyStr);
     //Ignoring third line
     getline(p3c_file,dummyStr);
     std::cout << dummyStr << std::endl;
     //go through N lines in the p3c file
     for (int i=0; i<N; i++) getline(p3c_file,dummyStr);
     
     //Ignoring first line
     getline(p3w_file, dummyStr); 
     std::cout << dummyStr << std::endl; 
     //Reading second line
     p3w_file >> dummy; set_t(dummy);
     p3w_file >> N;
     std::cout << "ignoring t= " << get_t() << ": loading " << N << " wall contacts ..." << std::endl;
     getline(p3w_file,dummyStr);
     //Ignoring third line
     getline(p3w_file,dummyStr);
     //go through N lines in the p3w file
     for (int i=0; i<N; i++) getline(p3w_file,dummyStr);
 
 }

template<StatType T>

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

protected

Definition at line 247 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);
     double nx, ny, nz;
     is >> dummy >> nx >> ny >> nz;
     double dx = (xmax-xmin)/(nx-1);
     double dy = (ymax-ymin)/(ny-1);
     double dz = (zmax-zmin)/(nz-1);
     VelocityProfile_D = Vec3D(dx,dy,dz);
     std::string statType;
     is >> dummy >> statType;
     getline(is,line_string);
     getline(is,line_string);
     
     set_nx(nx);
     set_ny(ny);
     set_nz(nz);
     set_xminStat(xmin);
     set_yminStat(ymin);
     set_zminStat(zmin);
     set_xmaxStat(xmax);
     set_ymaxStat(ymax);
     set_zmaxStat(zmax);
     
     VelocityProfile.resize(nx*ny*nz);
     Vec3D M;
     double rho;
     int n=0;
     for (int i=0; i<nx; i++)
     for (int j=0; j<ny; j++)
     for (int k=0; k<nz; k++) {
         is >> dummy >> dummy >> dummy >> dummy >> rho >> M.X >> M.Y >> M.Z;
         getline(is,line_string);
         VelocityProfile[n]=M/rho;
         n++;
     }
     
     //Close the file
     is.close();
     std::cout << "Loaded data file " << filename << std::endl;
     return true;
 }

template<StatType T>

void StatisticsVector< T >::output_statistics ( )

virtual

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

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

Todo:: {Particles shouldn't be unfixed, in case you do live statistics; please correct}

Todo:: We now have an idea of how to deal with fixed particles better, so this can be improved.

Reimplemented from MD.

Definition at line 1826 of file StatisticsVector.hcc.

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

template<StatType T>

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

Outputs member variable values to a std::string.

Outputs StatisticsVector.

Definition at line 584 of file StatisticsVector.hcc.

References Gaussian, HeavisideSphere, and Polynomial.

                                      {
     std::stringstream ss;
     ss  << "Statistical parameters: name=" << get_name() 
         << ",\n nx=" << nx
         << ", ny=" << ny
         << ", nz=" << nz 
         << ", CG_type=";
     if (get_CG_type()==HeavisideSphere) {
         ss << "HeavisideSphere";
     } else if (get_CG_type()==Gaussian) {
         ss << "Gaussian";
     } else if (get_CG_type()==Polynomial) {
         ss << get_PolynomialName();
     } else { std::cerr << "error in CG_function" << std::endl; exit(-1); }
 
     ss  << ", w=" << sqrt(w2)
         << ", cutoff=" << cutoff 
         << ",\n tStat=[" << get_tminStat() << "," << get_tmaxStat() << "]"
         << ", xStat=[" << get_xminStat() << "," << get_xmaxStat() << "]"
         << ", yStat=[" << get_yminStat() << "," << get_ymaxStat() << "]"
         << ", zStat=[" << get_zminStat() << "," << get_zmaxStat() << "]"
         << ",\n ignoreFixedParticles=" << ignoreFixedParticles
         << ", StressTypeForFixedParticles=" << StressTypeForFixedParticles
         << ",\n mirrorAtDomainBoundary=" << get_mirrorAtDomainBoundary()
         << ",\n doTimeAverage=" << get_doTimeAverage()
         << ", doVariance=" << get_doVariance()
         << ", doGradient=" << get_doGradient()
         << ", verbosity=" << verbosity
         << std::endl;
     return ss.str();
 }

template<StatType T>

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

Output coarse graining variables.

Output names of statistical variables.

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

Definition at line 619 of file StatisticsVector.hcc.

References Polynomial.

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

template<StatType T>

void StatisticsVector< T >::print_help ( )

Definition at line 476 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>

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

virtual

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

Reimplemented from MD.

Definition at line 736 of file StatisticsVector.hcc.

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

template<StatType T>

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

Definition at line 687 of file StatisticsVector.hcc.

References MD::read_next_from_data_file().

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

template<StatType T>

bool StatisticsVector< T >::read_next_from_p3p_file ( )

protected

Definition at line 968 of file StatisticsVector.hcc.

References constants::pi.

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

template<StatType T>

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

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

Definition at line 331 of file StatisticsVector.hcc.

References Gaussian, and HeavisideSphere.

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

template<StatType T>

void StatisticsVector< T >::reset_statistics ( )

virtual

Set all statistical variables to zero.

Definition at line 63 of file StatisticsVector.hcc.

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

template<StatType T>

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

protected

Definition at line 1511 of file StatisticsVector.hcc.

References BaseParticle::get_Position(), BaseParticle::get_Radius(), and Vec3D::Z.

                                                                     {
     return P->get_Radius()>rmin 
     && P->get_Radius()<rmax 
     && P->get_Position().Z<hmax;
 }

template<StatType T>

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

Definition at line 548 of file StatisticsVector.hcc.

References Gaussian, HeavisideSphere, and Polynomial.

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

template<StatType T>

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

Definition at line 577 of file StatisticsVector.hcc.

References Polynomial.

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

template<StatType T>

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

inline

Definition at line 143 of file StatisticsVector.h.

References StatisticsVector< T >::doGradient.

143 {doGradient = new_;};

StatisticsVector::doGradient

bool doGradient

Determines if gradient is outputted.

Definition: StatisticsVector.h:303

template<StatType T>

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

inline

Definition at line 128 of file StatisticsVector.h.

References StatisticsVector< T >::doTimeAverage.

128 {doTimeAverage = new_;};

StatisticsVector::doTimeAverage

bool doTimeAverage

Determines if output is averaged over time.

Definition: StatisticsVector.h:297

template<StatType T>

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

inline

Definition at line 232 of file StatisticsVector.h.

References StatisticsVector< T >::doublePoints.

232 {doublePoints=new_;}

StatisticsVector::doublePoints

bool doublePoints

Definition: StatisticsVector.h:363

template<StatType T>

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

inline

Definition at line 140 of file StatisticsVector.h.

References StatisticsVector< T >::doVariance.

140 {doVariance = new_;};

StatisticsVector::doVariance

bool doVariance

Determines if variance is outputted.

Definition: StatisticsVector.h:301

template<StatType T>

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

inline

Definition at line 81 of file StatisticsVector.h.

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

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

StatisticsVector::set_hx

void set_hx(Mdouble hx)

Definition: StatisticsVector.h:75

StatisticsVector::set_hy

void set_hy(Mdouble hy)

Definition: StatisticsVector.h:77

StatisticsVector::set_hz

void set_hz(Mdouble hz)

Definition: StatisticsVector.h:79

template<StatType T>

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

inline

Definition at line 242 of file StatisticsVector.h.

References StatisticsVector< T >::hmax.

242 {hmax=new_;}

StatisticsVector::hmax

Mdouble hmax

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

Definition: StatisticsVector.h:337

template<StatType T>

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

inline

Definition at line 75 of file StatisticsVector.h.

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

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

75 { set_nx( ceil((get_xmaxStat()-get_xminStat())/hx) ); };

StatisticsVector::get_xminStat

Mdouble get_xminStat()

Functions to acces and change the domain of statistics.

Definition: StatisticsVector.h:200

StatisticsVector::set_nx

void set_nx(int new_)

Definition: StatisticsVector.h:74

StatisticsVector::get_xmaxStat

Mdouble get_xmaxStat()

Definition: StatisticsVector.h:203

template<StatType T>

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

inline

Definition at line 77 of file StatisticsVector.h.

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

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

77 { set_ny( ceil((get_ymaxStat()-get_yminStat())/hy) ); };

StatisticsVector::get_yminStat

Mdouble get_yminStat()

Definition: StatisticsVector.h:201

StatisticsVector::set_ny

void set_ny(int new_)

Definition: StatisticsVector.h:76

StatisticsVector::get_ymaxStat

Mdouble get_ymaxStat()

Definition: StatisticsVector.h:204

template<StatType T>

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

inline

Definition at line 79 of file StatisticsVector.h.

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

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

79 { set_nz( ceil((get_zmaxStat()-get_zminStat())/hz) ); };

StatisticsVector::set_nz

void set_nz(int new_)

Definition: StatisticsVector.h:78

StatisticsVector::get_zmaxStat

Mdouble get_zmaxStat()

Definition: StatisticsVector.h:205

StatisticsVector::get_zminStat

Mdouble get_zminStat()

Definition: StatisticsVector.h:202

template<StatType T>

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

inline

Definition at line 149 of file StatisticsVector.h.

References StatisticsVector< T >::ignoreFixedParticles.

149 {ignoreFixedParticles = new_;};

StatisticsVector::ignoreFixedParticles

bool ignoreFixedParticles

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

Definition: StatisticsVector.h:344

template<StatType T>

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

inline

Definition at line 135 of file StatisticsVector.h.

References StatisticsVector< T >::StressTypeForFixedParticles.

135 {StressTypeForFixedParticles = 2+new_;};

StatisticsVector::StressTypeForFixedParticles

int StressTypeForFixedParticles

0 no Stress from fixed particles 1 Stress from fixed particles distributed between Contact and flowin...

Definition: StatisticsVector.h:349

template<StatType T>

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

inline

Definition at line 137 of file StatisticsVector.h.

References StatisticsVector< T >::mirrorAtDomainBoundary.

137 {mirrorAtDomainBoundary = new_;};

StatisticsVector::mirrorAtDomainBoundary

Mdouble mirrorAtDomainBoundary

Definition: StatisticsVector.h:356

template<StatType T>

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

inline

Definition at line 80 of file StatisticsVector.h.

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

80 { set_nx(n); set_ny(n); set_nz(n); };

StatisticsVector::set_nz

void set_nz(int new_)

Definition: StatisticsVector.h:78

StatisticsVector::set_nx

void set_nx(int new_)

Definition: StatisticsVector.h:74

StatisticsVector::set_ny

void set_ny(int new_)

Definition: StatisticsVector.h:76

template<StatType T>

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

inline

Definition at line 97 of file StatisticsVector.h.

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

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

StatisticsVector::nx

int nx

Grid size nx,ny,nz (by default the points of evaluation are placed in an grid on the domain [xminStat...

Definition: StatisticsVector.h:264

StatisticsVector::ny

int ny

see nx

Definition: StatisticsVector.h:266

StatisticsVector::nz

int nz

see nx

Definition: StatisticsVector.h:268

template<StatType T>

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

inline

Definition at line 74 of file StatisticsVector.h.

References StatisticsVector< T >::nx.

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

74 { nx = new_;};

StatisticsVector::nx

int nx

Grid size nx,ny,nz (by default the points of evaluation are placed in an grid on the domain [xminStat...

Definition: StatisticsVector.h:264

template<>

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

Definition at line 1955 of file StatisticsVector.hcc.

1955 { nx = 1; }

StatisticsVector::nx

int nx

Grid size nx,ny,nz (by default the points of evaluation are placed in an grid on the domain [xminStat...

Definition: StatisticsVector.h:264

template<>

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

Definition at line 1958 of file StatisticsVector.hcc.

1958 { nx = 1; }

StatisticsVector::nx

int nx

Grid size nx,ny,nz (by default the points of evaluation are placed in an grid on the domain [xminStat...

Definition: StatisticsVector.h:264

template<>

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

Definition at line 1960 of file StatisticsVector.hcc.

1960 { nx = 1; }

StatisticsVector::nx

int nx

Grid size nx,ny,nz (by default the points of evaluation are placed in an grid on the domain [xminStat...

Definition: StatisticsVector.h:264

template<>

void StatisticsVector< O >::set_nx ( int new_ UNUSED )

Definition at line 1962 of file StatisticsVector.hcc.

1962 { nx = 1; }

StatisticsVector::nx

int nx

Grid size nx,ny,nz (by default the points of evaluation are placed in an grid on the domain [xminStat...

Definition: StatisticsVector.h:264

template<StatType T>

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

inline

Definition at line 76 of file StatisticsVector.h.

References StatisticsVector< T >::ny.

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

76 { ny = new_;};

StatisticsVector::ny

int ny

see nx

Definition: StatisticsVector.h:266

template<>

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

Definition at line 1954 of file StatisticsVector.hcc.

1954 { ny = 1; }

StatisticsVector::ny

int ny

see nx

Definition: StatisticsVector.h:266

template<>

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

Definition at line 1956 of file StatisticsVector.hcc.

1956 { ny = 1; }

StatisticsVector::ny

int ny

see nx

Definition: StatisticsVector.h:266

template<>

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

Definition at line 1961 of file StatisticsVector.hcc.

1961 { ny = 1; }

StatisticsVector::ny

int ny

see nx

Definition: StatisticsVector.h:266

template<>

void StatisticsVector< O >::set_ny ( int new_ UNUSED )

Definition at line 1963 of file StatisticsVector.hcc.

1963 { ny = 1; }

StatisticsVector::ny

int ny

see nx

Definition: StatisticsVector.h:266

template<StatType T>

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

inline

Definition at line 78 of file StatisticsVector.h.

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

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

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

MD::get_dim

int get_dim()

Allows the dimension of the simulation to be accessed.

Definition: MD.h:369

StatisticsVector::nz

int nz

see nx

Definition: StatisticsVector.h:268

template<>

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

Definition at line 1953 of file StatisticsVector.hcc.

1953 { nz = 1; }

StatisticsVector::nz

int nz

see nx

Definition: StatisticsVector.h:268

template<>

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

Definition at line 1957 of file StatisticsVector.hcc.

1957 { nz = 1; }

StatisticsVector::nz

int nz

see nx

Definition: StatisticsVector.h:268

template<>

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

Definition at line 1959 of file StatisticsVector.hcc.

1959 { nz = 1; }

StatisticsVector::nz

int nz

see nx

Definition: StatisticsVector.h:268

template<>

void StatisticsVector< O >::set_nz ( int new_ UNUSED )

Definition at line 1964 of file StatisticsVector.hcc.

1964 { nz = 1; }

StatisticsVector::nz

int nz

see nx

Definition: StatisticsVector.h:268

template<StatType T>

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

inline

Definition at line 159 of file StatisticsVector.h.

References StatisticsVector< T >::periodicWalls.

Referenced by Statistics().

159 {periodicWalls = new_;};

StatisticsVector::periodicWalls

bool periodicWalls

Turns off periodic walls before evaluation (needed for averaging, because we do not yet check if part...

Definition: StatisticsVector.h:342

template<StatType T>

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

inline

Definition at line 216 of file StatisticsVector.h.

References StatisticsVector< T >::CGPolynomial.

                                                                                    {    
         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 220 of file StatisticsVector.h.

References StatisticsVector< T >::CGPolynomial.

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

template<StatType T>

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

inline

Definition at line 224 of file StatisticsVector.h.

References StatisticsVector< T >::CGPolynomial.

                                                   { 
         CGPolynomial.set_name(new_name);
     }

template<StatType T>

void StatisticsVector< T >::set_Positions ( )

Set position of StatisticsPoint points and set variables to 0.

Set position of statistics points and set variables to 0.

Definition at line 1365 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 >::initialize_statistics().

 {
     //Dimensions of the domain / n 
     Vec3D diff = Vec3D(
         (get_xmaxStat()-get_xminStat())/nx, 
         (get_ymaxStat()-get_yminStat())/ny,
         (get_zmaxStat()-get_zminStat())/nz);
     //Center of the domain
     Vec3D avg = 0.5 * Vec3D(
         get_xmaxStat()+get_xminStat(), 
         get_ymaxStat()+get_yminStat(), 
         get_zmaxStat()+get_zminStat());
     //Left endpoint of the domain
     Vec3D Min = Vec3D(
         get_xminStat(), 
         get_yminStat(), 
         get_zminStat());
 
     int num[3]={0,0,0};
     if (get_mirrorAtDomainBoundary()) {
         if (statType==X||statType==XY||statType==XZ||statType==XYZ) {
             //add so many points to domain on both ends
             //todo: the 8 is arbitrary
             num[0] = floor(std::max(get_mirrorAtDomainBoundary(),get_cutoff())/diff.X);
             nx+=2*num[0];
             Min.X-=num[0]*diff.X;
         }
         if (statType==Y||statType==XY||statType==YZ||statType==XYZ) {
             num[1] = floor(std::max(get_mirrorAtDomainBoundary(),get_cutoff())/diff.Y);
             ny+=2*num[1];
             Min.Y-=num[1]*diff.Y;
         }
         if (statType==Z||statType==XZ||statType==XZ||statType==XYZ) {
             num[2] = floor(std::max(get_mirrorAtDomainBoundary(),get_cutoff())/diff.Z);
             nz+=2*num[2];
             Min.Z-=num[2]*diff.Z;
         }
     }
 
     //Set position of statistics points and set variables to 0
     int N =std::max(nx,1)*std::max(ny,1)*std::max(nz,1);
     Points.resize(N);
     int n = 0;
     for (int i=0; i<std::max(nx,1); i++)
     for (int j=0; j<std::max(ny,1); j++)
     for (int k=0; k<std::max(nz,1); k++) {
         Points[n].set_Position( Vec3D( (nx>1)?(Min.X+diff.X*(i+0.5)):avg.X, 
             (ny>1)?(Min.Y+diff.Y*(j+0.5)):avg.Y, 
             (nz>1)?(Min.Z+diff.Z*(k+0.5)):avg.Z ));
         if (doublePoints) Points[n].set_Position( Points[n].get_Position()
             -Vec3D((i%2)?(0.99*diff.X):0,(j%2)?(0.99*diff.Y):0,(k%2)?(0.99*diff.Z):0) );
         Points[n].set_CG_invvolume();
         Points[n].set_zero();
         if (get_mirrorAtDomainBoundary())
         if ((i<num[0])||(i>nx-num[0])||(j<num[1])||(j>ny-num[1])||(k<num[2])||(k>nz-num[2])) { 
             Points[n].mirrorParticle = n + 
                 + ((i<num[0])?(2*(num[0]-i)-1)*ny*nz:0)+((i>nx-num[0])?(2*(nx-num[0]-i)+1)*ny*nz:0)
                 + ((j<num[1])?(2*(num[1]-j)-1)*nz:0   )+((j>ny-num[1])?(2*(ny-num[1]-j)+1)*nz:0)
                 + ((k<num[2])?(2*(num[2]-k)-1):0      )+((k>nz-num[2])?(2*(nz-num[2]-k)+1):0); 
         }
         n++;
     }
     if (statType==RAZ || statType==RZ  || statType==RA  || statType==AZ || statType==R || statType==A) {
         for (unsigned int k=0; k<Points.size(); k++) {
             Points[k].set_Position(Points[k].get_Position().GetFromCylindricalCoordinates());
         }
     }
     if (get_doGradient()) {
         dx.resize(N);
         dy.resize(N);
         dz.resize(N);
         for (int n=0; n<N; n++) {
             dx[n].set_Position(Points[n].get_Position());
             dx[n].mirrorParticle=Points[n].mirrorParticle;
             dx[n].set_zero();
             dy[n].set_Position(Points[n].get_Position());
             dy[n].mirrorParticle=Points[n].mirrorParticle;
             dy[n].set_zero();
             dz[n].set_Position(Points[n].get_Position());
             dz[n].mirrorParticle=Points[n].mirrorParticle;
             dz[n].set_zero();
         }
     }
 
     //Set TimeAverage the same way
     if (get_doTimeAverage()) {
         timeAverage.resize(N);
         int n = 0;
         for (int i=0; i<std::max(nx,1); i++)
             for (int j=0; j<std::max(ny,1); j++)
                 for (int k=0; k<std::max(nz,1); k++) {
                     timeAverage[n].set_Position(Points[n].get_Position());
                     timeAverage[n].set_zero();
                     timeAverage[n].mirrorParticle=Points[n].mirrorParticle;
                     n++;
                 }
 
         if (get_doVariance()) {
             timeVariance.resize(N);
             for (int n=0; n<N; n++) {
                 timeVariance[n].set_Position(timeAverage[n].get_Position());
                 timeVariance[n].mirrorParticle=Points[n].mirrorParticle;
                 timeVariance[n].set_zero();
             }
         }
 
         if (get_doGradient()) {
             dxTimeAverage.resize(N);
             dyTimeAverage.resize(N);
             dzTimeAverage.resize(N);
             for (int n=0; n<N; n++) {
                 dxTimeAverage[n].set_Position(timeAverage[n].get_Position());
                 dxTimeAverage[n].mirrorParticle=Points[n].mirrorParticle;
                 dxTimeAverage[n].set_zero();
                 dyTimeAverage[n].set_Position(timeAverage[n].get_Position());
                 dyTimeAverage[n].mirrorParticle=Points[n].mirrorParticle;
                 dyTimeAverage[n].set_zero();
                 dzTimeAverage[n].set_Position(timeAverage[n].get_Position());
                 dzTimeAverage[n].mirrorParticle=Points[n].mirrorParticle;
                 dzTimeAverage[n].set_zero();
             }
         }
 
     } //end if (get_doTimeAverage()) 
 }

template<StatType T>

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

inline

Definition at line 241 of file StatisticsVector.h.

References StatisticsVector< T >::rmax.

241 {rmax=new_;}

StatisticsVector::rmax

Mdouble rmax

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

Definition: StatisticsVector.h:335

template<StatType T>

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

inline

Definition at line 240 of file StatisticsVector.h.

References StatisticsVector< T >::rmin.

240 {rmin=new_;}

StatisticsVector::rmin

Mdouble rmin

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

Definition: StatisticsVector.h:333

template<StatType T>

void StatisticsVector< T >::set_statType ( )

template<>

void StatisticsVector< RAZ >::set_statType ( )

Definition at line 1966 of file StatisticsVector.hcc.

References RAZ.

1966 { statType=RAZ; }

RAZ

Definition: StatisticsVector.h:34

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

template<>

void StatisticsVector< RZ >::set_statType ( )

Definition at line 1967 of file StatisticsVector.hcc.

References RZ.

1967 { statType=RZ; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

RZ

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< AZ >::set_statType ( )

Definition at line 1968 of file StatisticsVector.hcc.

References AZ.

1968 { statType=AZ; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

AZ

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< RA >::set_statType ( )

Definition at line 1969 of file StatisticsVector.hcc.

References RA.

1969 { statType=RA; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

RA

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< A >::set_statType ( )

Definition at line 1970 of file StatisticsVector.hcc.

References A.

1970 { statType=A; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

A

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< R >::set_statType ( )

Definition at line 1971 of file StatisticsVector.hcc.

References R.

1971 { statType=R; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

R

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< XYZ >::set_statType ( )

Definition at line 1972 of file StatisticsVector.hcc.

References XYZ.

1972 { statType=XYZ; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

XYZ

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< XY >::set_statType ( )

Definition at line 1973 of file StatisticsVector.hcc.

References XY.

1973 { statType=XY; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

XY

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< XZ >::set_statType ( )

Definition at line 1974 of file StatisticsVector.hcc.

References XZ.

1974 { statType=XZ; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

XZ

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< YZ >::set_statType ( )

Definition at line 1975 of file StatisticsVector.hcc.

References YZ.

1975 { statType=YZ; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

YZ

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< X >::set_statType ( )

Definition at line 1976 of file StatisticsVector.hcc.

References X.

1976 { statType=X; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

X

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< Y >::set_statType ( )

Definition at line 1977 of file StatisticsVector.hcc.

References Y.

1977 { statType=Y; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

Y

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< Z >::set_statType ( )

Definition at line 1978 of file StatisticsVector.hcc.

References Z.

1978 { statType=Z; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

Z

Definition: StatisticsVector.h:34

template<>

void StatisticsVector< O >::set_statType ( )

Definition at line 1979 of file StatisticsVector.hcc.

References O.

1979 { statType=O; }

StatisticsVector::statType

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

Definition: StatisticsVector.h:262

O

Definition: StatisticsVector.h:34

template<StatType T>

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

inline

Definition at line 131 of file StatisticsVector.h.

References StatisticsVector< T >::StressTypeForFixedParticles.

131 {StressTypeForFixedParticles = new_;};

StatisticsVector::StressTypeForFixedParticles

int StressTypeForFixedParticles

0 no Stress from fixed particles 1 Stress from fixed particles distributed between Contact and flowin...

Definition: StatisticsVector.h:349

template<StatType T>

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

inline

Definition at line 146 of file StatisticsVector.h.

References StatisticsVector< T >::superexact.

146 {superexact = new_;};

StatisticsVector::superexact

bool superexact

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

Definition: StatisticsVector.h:360

template<StatType T>

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

inline

Definition at line 236 of file StatisticsVector.h.

References StatisticsVector< T >::nTimeAverageReset.

236 {nTimeAverageReset=new_;}

StatisticsVector::nTimeAverageReset

int nTimeAverageReset

Determines after how many timesteps the time average is reset.

Definition: StatisticsVector.h:299

template<StatType T>

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

inline

Definition at line 87 of file StatisticsVector.h.

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

87 {tintStat=t;};

StatisticsVector::tintStat

Mdouble tintStat

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

Definition: StatisticsVector.h:328

MD::t

Mdouble t

This stores the current time.

Definition: MD.h:687

template<StatType T>

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

inline

Definition at line 86 of file StatisticsVector.h.

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

86 {tmaxStat=t;};

StatisticsVector::tmaxStat

Mdouble tmaxStat

Statistical output will only be created if t

Definition: StatisticsVector.h:326

MD::t

Mdouble t

This stores the current time.

Definition: MD.h:687

template<StatType T>

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

inline

Definition at line 85 of file StatisticsVector.h.

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

85 {tminStat=t;};

StatisticsVector::tminStat

Mdouble tminStat

Statistical output will only be created if t>tminStat.

Definition: StatisticsVector.h:324

MD::t

Mdouble t

This stores the current time.

Definition: MD.h:687

template<StatType T>

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

inline

Definition at line 153 of file StatisticsVector.h.

References StatisticsVector< T >::verbosity.

153 {verbosity = new_;};

StatisticsVector::verbosity

int verbosity

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

Definition: StatisticsVector.h:354

template<StatType T>

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

inline

Set CG variables w2 and CG_invvolume.

Definition at line 101 of file StatisticsVector.h.

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

101 {set_w2(sqr(w));};

sqr

#define sqr(a)

Definition: ExtendedMath.h:36

StatisticsVector::set_w2

void set_w2(Mdouble new_)

Set CG variables w2 and CG_invvolume.

Definition: StatisticsVector.hcc:46

template<StatType T>

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

Set CG variables w2 and CG_invvolume.

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

Definition at line 46 of file StatisticsVector.hcc.

References Gaussian, HeavisideSphere, Polynomial, and sqr.

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

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

template<StatType T>

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

inline

Definition at line 162 of file StatisticsVector.h.

References StatisticsVector< T >::w_over_rmax.

162 {w_over_rmax = new_;};

StatisticsVector::w_over_rmax

Mdouble w_over_rmax

if w is not set manually then w will be set by multiplying this value by the largest particle radius ...

Definition: StatisticsVector.h:319

template<StatType T>

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

inline

Definition at line 156 of file StatisticsVector.h.

References StatisticsVector< T >::walls.

156 {walls = new_;};

StatisticsVector::walls

bool walls

Turns off walls before evaluation.

Definition: StatisticsVector.h:339

template<StatType T>

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

inline

Definition at line 209 of file StatisticsVector.h.

References StatisticsVector< T >::xmaxStat.

209 {xmaxStat=xmaxStat_;};

StatisticsVector::xmaxStat

Mdouble xmaxStat

Definition: StatisticsVector.h:271

template<StatType T>

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

inline

Definition at line 206 of file StatisticsVector.h.

References StatisticsVector< T >::xminStat.

206 {xminStat=xminStat_;};

StatisticsVector::xminStat

Mdouble xminStat

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

Definition: StatisticsVector.h:271

template<StatType T>

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

inline

Definition at line 210 of file StatisticsVector.h.

References StatisticsVector< T >::ymaxStat.

210 {ymaxStat=ymaxStat_;};

StatisticsVector::ymaxStat

Mdouble ymaxStat

Definition: StatisticsVector.h:271

template<StatType T>

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

inline

Definition at line 207 of file StatisticsVector.h.

References StatisticsVector< T >::yminStat.

207 {yminStat=yminStat_;};

StatisticsVector::yminStat

Mdouble yminStat

Definition: StatisticsVector.h:271

template<StatType T>

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

inline

Definition at line 211 of file StatisticsVector.h.

References StatisticsVector< T >::zmaxStat.

211 {zmaxStat=zmaxStat_;};

StatisticsVector::zmaxStat

Mdouble zmaxStat

Definition: StatisticsVector.h:271

template<StatType T>

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

inline

Definition at line 208 of file StatisticsVector.h.

References StatisticsVector< T >::zminStat.

208 {zminStat=zminStat_;};

StatisticsVector::zminStat

Mdouble zminStat

Definition: StatisticsVector.h:271

template<StatType T>

Mdouble StatisticsVector< T >::setInfinitelyLongDistance ( )

template<>

Mdouble StatisticsVector< RAZ >::setInfinitelyLongDistance ( )

todo{check}

Definition at line 1981 of file StatisticsVector.hcc.

                                                                     {
     return std::min(std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
                    std::max((get_xmaxStat()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X)),
                    std::max((get_ymaxStat()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y));
 }

template<>

Mdouble StatisticsVector< RA >::setInfinitelyLongDistance ( )

Definition at line 1987 of file StatisticsVector.hcc.

                                                                    {
     return std::min(std::max((get_ymaxStat()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y),
                std::max((get_xmaxStat()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X));
 }

template<>

Mdouble StatisticsVector< RZ >::setInfinitelyLongDistance ( )

Definition at line 1991 of file StatisticsVector.hcc.

                                                                    {
     return std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
                std::max((get_xmaxStat()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X));
 }

template<>

Mdouble StatisticsVector< AZ >::setInfinitelyLongDistance ( )

Definition at line 1995 of file StatisticsVector.hcc.

                                                                    {
     return std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
                std::max((get_ymaxStat()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y));
 }

template<>

Mdouble StatisticsVector< R >::setInfinitelyLongDistance ( )

Definition at line 1999 of file StatisticsVector.hcc.

                                                                   {
     return std::max((get_xmaxStat()-P2.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P2.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X);
 }

template<>

Mdouble StatisticsVector< A >::setInfinitelyLongDistance ( )

Definition at line 2002 of file StatisticsVector.hcc.

                                                                   {
     return std::max((get_ymaxStat()-P2.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P2.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y);
 }

template<>

Mdouble StatisticsVector< XYZ >::setInfinitelyLongDistance ( )

Definition at line 2007 of file StatisticsVector.hcc.

                                                                     {
     return std::min(std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
                    std::max((get_xmaxStat()-P1.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P1.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X)),
                    std::max((get_ymaxStat()-P1.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P1.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y));
 }

template<>

Mdouble StatisticsVector< XY >::setInfinitelyLongDistance ( )

Definition at line 2012 of file StatisticsVector.hcc.

                                                                    {
     return std::min(std::max((get_ymaxStat()-P1.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P1.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y),
                std::max((get_xmaxStat()-P1.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P1.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X));
 }

template<>

Mdouble StatisticsVector< XZ >::setInfinitelyLongDistance ( )

Definition at line 2016 of file StatisticsVector.hcc.

                                                                    {
     return std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
                std::max((get_xmaxStat()-P1.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P1.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X));
 }

template<>

Mdouble StatisticsVector< YZ >::setInfinitelyLongDistance ( )

Definition at line 2020 of file StatisticsVector.hcc.

                                                                    {
     return std::min(fabs((P1.Z-get_zminStat()+5*get_w())/P1_P2_normal.Z),
                std::max((get_ymaxStat()-P1.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P1.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y));
 }

template<>

Mdouble StatisticsVector< X >::setInfinitelyLongDistance ( )

Definition at line 2024 of file StatisticsVector.hcc.

                                                                   {
     return std::max((get_xmaxStat()-P1.X+5*sqrt(get_w2()))/P1_P2_normal.X,-(P1.X-get_xminStat()+5*sqrt(get_w2()))/P1_P2_normal.X);
 }

template<>

Mdouble StatisticsVector< Y >::setInfinitelyLongDistance ( )

Definition at line 2027 of file StatisticsVector.hcc.

                                                                   {
     return std::max((get_ymaxStat()-P1.Y+5*sqrt(get_w2()))/P1_P2_normal.Y,-(P1.Y-get_yminStat()+5*sqrt(get_w2()))/P1_P2_normal.Y);
 }

template<>

Mdouble StatisticsVector< Z >::setInfinitelyLongDistance ( )

Definition at line 2030 of file StatisticsVector.hcc.

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

template<>

double StatisticsVector< O >::setInfinitelyLongDistance ( )

Definition at line 2033 of file StatisticsVector.hcc.

                                                                  {
     return 0;
 }

template<StatType T>

void StatisticsVector< T >::statistics_from_fstat_and_data ( )

get StatisticsPoint

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

Definition at line 860 of file StatisticsVector.hcc.

References MD::print().

Referenced by Statistics().

 {
     //This opens the file the data will be recalled from
     set_data_filename();
     open_data_file(std::fstream::in);
 
     //load first set of timesteps to obtain xdim, ydim, zdim
 
     //try to read data files; if data.0000 does not exist, try the next, up to .9999
     for (int i=1; i<10000; i++) {
         if(read_next_from_data_file(format)) break;
         set_file_counter(i);
     }
 
     //set tminstat to smallest time evaluated
     if (get_t()>=get_tminStat()) set_tminStat(get_t()-get_dt());
                 
     //This opens the file the force data will be recalled from
     set_fstat_filename();
     open_fstat_file(std::fstream::in);
 
     //This creates the file statistics will be saved to
     //open_stat_file(std::fstream::out);
 
     if (verbosity>1) {
         std::cout << "Input from " << get_data_filename() << std::endl;
         std::cout << "Input from " << get_fstat_filename() << std::endl;
         std::cout << "Output to " << get_stat_filename() << std::endl << std::endl;
     }
     
     initialize_statistics();
     
     //std::couts collision time
     Mdouble mass = get_Mass_from_Radius(getLargestParticle()->get_Radius());
     if (verbosity>1) std::cout << "Collision Time " << get_collision_time(mass) 
         << " and restitution coefficient " << get_restitution_coefficient(mass) 
         << " for mass " << mass << std::endl
         << std::endl;
     //std::couts particles
     if (verbosity>2) {MD::print(std::cout,true); std::cout << std::endl;}
     else if (verbosity) {MD::print(std::cout); std::cout << std::endl;}
     
     //this increases the file counter
     if (get_options_data()==2 && get_t()<get_tminStat()) {
         if (verbosity>1) find_next_data_file(get_tminStat());
         else find_next_data_file(get_tminStat(), false);
         for(int i=1;i<get_file_counter();i++) jump_fstat();
         
     }
     
     //Output statistics for each time step 
     std::cout<<"Start statistics"<<std::endl;
     do {
         if (get_t()>get_tmaxStat()) { 
             std::cout << "reached end t="<<std::setprecision (9) <<get_t()<<" tmaxStat=" << std::setprecision (9) <<get_tmaxStat()<< std::endl; break;
         } else if (get_t()>=get_tminStat()) {
             if (verbosity>1) std::cout << "Get statistics for t=" << std::setprecision (9) << get_t() << std::endl;
             else if (verbosity) std::cout << "Get statistics for t=" << std::setprecision (9) << get_t() << " \r";
             if (verbosity>1) std::cout << "#particles="<< getParticleHandler().getNumberOfObjects() << std::endl;
             
             if (!walls) getWallHandler().clear();
             if (!periodicWalls) getBoundaryHandler().clear();
                     
             output_statistics();
             gather_force_statistics_from_fstat_and_data();
             process_statistics(true);
 
         } else {
             if (verbosity>1) std::cout<<"Jumping statistics t="<<get_t()<<" tminStat()="<<get_tminStat()<<std::endl;
             jump_fstat();
         }
         //if step size>1, skip a number of steps.
         if (get_options_data()!=2) for (unsigned int i=1; i<get_step_size(); i++) {
             //you only get here if you change the step size and you use a single data file
             read_next_from_data_file(format);
             jump_fstat();
         }
     } while (read_next_from_data_file(format));
     //set tmax to largest time evaluated
     if (get_t()<get_tmaxStat()) set_tmaxStat(get_t());
 
     finish_statistics();
     get_data_file().close();
     get_fstat_file().close();
 }

template<StatType T>

void StatisticsVector< T >::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 1033 of file StatisticsVector.hcc.

References MD::create_xballs_script(), MD::output_xballs_data(), MD::print(), MD::save_restart_data(), MD::set_dim(), and MD::set_format().

 {
     //This opens the files the data will be recalled from
     std::string p3p_filename = problem_name.str().c_str(); p3p_filename += ".p3p";
     bool opened = open_file (p3p_file, p3p_filename, 1,std::fstream::in);
     std::string p3c_filename = problem_name.str().c_str(); p3c_filename += ".p3c";
     open_file (p3c_file, p3c_filename, 1,std::fstream::in);
     std::string p3w_filename = problem_name.str().c_str(); p3w_filename += ".p3w";
     open_file (p3w_file, p3w_filename, 1,std::fstream::in);
     int version=3;
     //If p3 is not available, open p4 files
     if (!opened) {
         std::cerr << "Warning could not open " << p3p_filename << std::endl; 
         version=4;
         p3p_filename = problem_name.str().c_str(); p3p_filename += ".p4p";
         opened = open_file (p3p_file, p3p_filename, 1,std::fstream::in);
         p3c_filename = problem_name.str().c_str(); p3c_filename += ".p4c";
         opened = open_file (p3c_file, p3c_filename, 1,std::fstream::in);
         p3w_filename = problem_name.str().c_str(); p3w_filename += ".p4w";
         opened = open_file (p3w_file, p3w_filename, 1,std::fstream::in);
         if (!opened) {
             std::cerr << "could not open " << p3p_filename << std::endl; 
             exit(-1);
         }
     }
     std::cout << "opened " << p3p_filename << std::endl;
 
     //load first set of timesteps to obtain xdim, ydim, zdim
 
     //try to read data files; if data.0000 does not exist, try the next, up to .9999
     read_next_from_p3p_file();
 
     //set tminstat to smallest time evaluated
     if (get_t()>=get_tminStat()) set_tminStat(get_t()-get_dt());
     MD::save_restart_data();
 
     set_data_filename();
     open_data_file(std::fstream::out);
     MD::set_format(14);
     MD::set_dim(3);
     Species[0].set_dim_particle(3);
     MD::output_xballs_data();
     MD::create_xballs_script();
     //exit(-1);
                 
     //This creates the file statistics will be saved to
     //open_stat_file(std::fstream::out);
 
     initialize_statistics();
     
     //std::couts collision time
     Mdouble mass = get_Mass_from_Radius(getLargestParticle()->get_Radius());
     if (verbosity>1) std::cout << "Collision Time " << get_collision_time(mass) 
         << " and restitution coefficient " << get_restitution_coefficient(mass) 
         << " for mass " << mass << std::endl
         << std::endl;
     //std::couts particles
     if (verbosity>2) {MD::print(std::cout,true); std::cout << std::endl;}
     else if (verbosity) {MD::print(std::cout); std::cout << std::endl;}
         
     //Output statistics for each time step 
     std::cout<<"Start statistics"<<std::endl;
     do {
         if (get_t()>get_tmaxStat()) { 
             std::cout << "reached end t="<<std::setprecision (9) <<get_t() <<" tmaxStat=" << std::setprecision (9) <<get_tmaxStat()<< std::endl; 
             break;
         } else if (get_t()>=get_tminStat()) {
             if (verbosity>1) std::cout << "Get statistics for t=" << std::setprecision (9) << get_t() << std::endl;
             else if (verbosity) std::cout << "Get statistics for t=" << std::setprecision (9) << get_t() << " \r";
             if (verbosity>1) std::cout << "#particles="<< getParticleHandler().getNumberOfObjects() << std::endl;
             
             if (!walls) getWallHandler().clear();
             if (!periodicWalls) getBoundaryHandler().clear();
                     
             output_statistics();
             gather_force_statistics_from_p3c(version);
             process_statistics(true);
 
         } else {
             if (verbosity>1) std::cout 
             << "Jumping statistics t=" << get_t() 
             << " tminStat()=" << get_tminStat() << std::endl;
             jump_p3c();
         }
         //if step size>1, skip a number of steps.
         for (unsigned int i=1; i<get_step_size(); i++) {
             //you only get here if you change the step size and you use a single data file
             read_next_from_p3p_file();
             jump_p3c();
         }
     } while (read_next_from_p3p_file());
     //set tmax to largest time evaluated
     if (get_t()<get_tmaxStat()) set_tmaxStat(get_t());
 
     finish_statistics();
     p3p_file.close();
     p3c_file.close();
     p3w_file.close();
 }

template<StatType T>

void StatisticsVector< T >::verbose ( )

inline

Definition at line 152 of file StatisticsVector.h.

References StatisticsVector< T >::verbosity.

152 {verbosity = 2;};

StatisticsVector::verbosity

int verbosity

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

Definition: StatisticsVector.h:354

template<StatType T>

void StatisticsVector< T >::write_statistics ( )

Writes regular statistics.

Definition at line 705 of file StatisticsVector.hcc.

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

template<StatType T>

void StatisticsVector< T >::write_time_average_statistics ( )

Writes out time averaged statistics.

Definition at line 788 of file StatisticsVector.hcc.

References StatisticsPoint< T >::print().

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

Member Data Documentation

template<StatType T>

CG StatisticsVector< T >::CG_type

protected

coarse graining type (Gaussian, Heaviside, Polynomial)

Definition at line 309 of file StatisticsVector.h.

Referenced by StatisticsVector< T >::get_CG_type(), 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 311 of file StatisticsVector.h.

Referenced by StatisticsVector< T >::evaluateIntegral(), StatisticsVector< T >::evaluatePolynomial(), StatisticsVector< T >::evaluatePolynomialGradient(), StatisticsVector< T >::get_PolynomialName(), StatisticsVector< T >::set_Polynomial(), StatisticsVector< T >::set_PolynomialName(), 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 317 of file StatisticsVector.h.

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

template<StatType T>

Mdouble StatisticsVector< T >::cutoff2

protected

Definition at line 317 of file StatisticsVector.h.

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

template<StatType T>

bool StatisticsVector< T >::doGradient

protected

Determines if gradient is outputted.

Definition at line 303 of file StatisticsVector.h.

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

template<StatType T>

bool StatisticsVector< T >::doTimeAverage

protected

Determines if output is averaged over time.

Definition at line 297 of file StatisticsVector.h.

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

template<StatType T>

bool StatisticsVector< T >::doublePoints

protected

Definition at line 363 of file StatisticsVector.h.

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

template<StatType T>

bool StatisticsVector< T >::doVariance

protected

Determines if variance is outputted.

Definition at line 301 of file StatisticsVector.h.

Referenced by StatisticsVector< T >::get_doVariance(), StatisticsVector< T >::set_doVariance(), 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 279 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 291 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 281 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 293 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 283 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 295 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::format

protected

Definition at line 355 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 337 of file StatisticsVector.h.

Referenced by StatisticsVector< T >::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 344 of file StatisticsVector.h.

Referenced by StatisticsVector< T >::get_ignoreFixedParticles(), StatisticsVector< T >::set_ignoreFixedParticles(), 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 330 of file StatisticsVector.h.

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

template<StatType T>

bool StatisticsVector< T >::isMDCLR

protected

Definition at line 358 of file StatisticsVector.h.

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

template<StatType T>

Mdouble StatisticsVector< T >::mirrorAtDomainBoundary

protected

Definition at line 356 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::nTimeAverage

protected

A counter needed to average over time.

Definition at line 305 of file StatisticsVector.h.

Referenced by StatisticsVector< T >::get_nTimeAverage(), 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 299 of file StatisticsVector.h.

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

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 264 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::nxMirrored

protected

extension of grid size from mirrored points

Definition at line 273 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::ny

protected

see nx

Definition at line 266 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::nyMirrored

protected

Definition at line 273 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::nz

protected

see nx

Definition at line 268 of file StatisticsVector.h.

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

template<StatType T>

int StatisticsVector< T >::nzMirrored

protected

Definition at line 273 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 370 of file StatisticsVector.h.

template<StatType T>

Vec3D StatisticsVector< T >::P1_P2_CollisionalHeatFlux

protected

not yet working

Definition at line 391 of file StatisticsVector.h.

template<StatType T>

Vec3D StatisticsVector< T >::P1_P2_Contact

protected

Definition at line 381 of file StatisticsVector.h.

template<StatType T>

Matrix3D StatisticsVector< T >::P1_P2_ContactCoupleStress

protected

Definition at line 380 of file StatisticsVector.h.

template<StatType T>

Mdouble StatisticsVector< T >::P1_P2_Dissipation

protected

not yet working

Definition at line 393 of file StatisticsVector.h.

template<StatType T>

Mdouble StatisticsVector< T >::P1_P2_distance

protected

Length of contact line.

Definition at line 376 of file StatisticsVector.h.

template<StatType T>

MatrixSymmetric3D StatisticsVector< T >::P1_P2_Fabric

protected

Fabric.

Definition at line 389 of file StatisticsVector.h.

template<StatType T>

Vec3D StatisticsVector< T >::P1_P2_normal

protected

Direction of contact.

Definition at line 374 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 378 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 385 of file StatisticsVector.h.

template<StatType T>

Mdouble StatisticsVector< T >::P1_P2_Potential

protected

not yet working

Definition at line 395 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 383 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 387 of file StatisticsVector.h.

template<StatType T>

Vec3D StatisticsVector< T >::P2

protected

Position of second contact point.

Definition at line 372 of file StatisticsVector.h.

template<StatType T>

std::fstream StatisticsVector< T >::p3c_file

protected

Definition at line 405 of file StatisticsVector.h.

template<StatType T>

std::fstream StatisticsVector< T >::p3p_file

protected

Definition at line 404 of file StatisticsVector.h.

template<StatType T>

std::fstream StatisticsVector< T >::p3w_file

protected

Definition at line 406 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 342 of file StatisticsVector.h.

Referenced by StatisticsVector< T >::get_periodicWalls(), StatisticsVector< T >::set_periodicWalls(), 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 277 of file StatisticsVector.h.

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

template<StatType T>

Mdouble StatisticsVector< T >::rmax

protected

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

Definition at line 335 of file StatisticsVector.h.

Referenced by StatisticsVector< T >::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 333 of file StatisticsVector.h.

Referenced by StatisticsVector< T >::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 262 of file StatisticsVector.h.

template<StatType T>

int StatisticsVector< T >::StressTypeForFixedParticles

protected

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 349 of file StatisticsVector.h.

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

template<StatType T>

bool StatisticsVector< T >::superexact

protected

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

Definition at line 360 of file StatisticsVector.h.

Referenced by StatisticsVector< T >::get_superexact(), StatisticsVector< T >::set_superexact(), 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.