StatisticsVector.h

Go to the documentation of this file.

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#ifndef STATISTICSVECTOR_H
#define STATISTICSVECTOR_H
 
#include "DPMBase.h"
#include "Math/Matrix.h"
#include "Math/MatrixSymmetric.h"
#include <string.h>
#include <fstream>
#include <math.h>
#include "Species/LinearViscoelasticSpecies.h"
 
enum StatType
{
    O, X, Y, Z, XY, XZ, YZ, XYZ, RAZ, RZ, AZ, RA, R, A
};
 
#include "Math/NormalisedPolynomial.h"
#include "StatisticsPoint.h"
 
template<StatType T>
class StatisticsVector : public virtual DPMBase
{
public:
    void constructor();
    
    void constructor(std::string name);
    
    StatisticsVector()
    {
        constructor();
    }
    
    StatisticsVector(std::string name)
    {
        constructor(name);
    }
    
    StatisticsVector(StatisticsVector& other);
    
    StatisticsVector(int argc, char* argv[]);
    
    void readStatArguments(int argc, char* argv[]);
    
    std::string printStat();
    
    void set_statType();
    
    void writeOutputFiles();
    
    void print_help();
    
    //set functions should be called before executing the statistical routine 
    
    void setNX(int new_)
    {
        nx = new_;
    }
    
    void set_hx(Mdouble hx)
    {
        setNX(static_cast<int>(std::ceil((getXMaxStat() - getXMinStat()) / hx)));
    }
    
    void setNY(int new_)
    {
        ny = new_;
    }
    
    void set_hy(Mdouble hy)
    {
        setNY(static_cast<int>(std::ceil((getYMaxStat() - getYMinStat()) / hy)));
    }
    
    void setNZ(int new_)
    {
        nz = new_;
        if (getSystemDimensions() < 3)
        {
            logger(WARN, "Warning in set_nz: dimension less than 3");
        }
    }
    
    void set_hz(Mdouble hz)
    {
        setNZ(static_cast<int>(std::ceil((getZMaxStat() - getZMinStat()) / hz)));
    }
    
    void setN(int n)
    {
        setNX(n);
        setNY(n);
        setNZ(n);
    }
    
    void set_h(Mdouble h)
    {
        set_hx(h);
        set_hy(h);
        set_hz(h);
    }
    
    int getNX()
    {
        return nx;
    }
    
    int getNY()
    {
        return ny;
    }
    
    int getNZ()
    {
        return nz;
    }
    
    void setCGTimeMin(Mdouble t)
    {
        tMinStat = t;
    }
    
    void setTimeMaxStat(Mdouble t)
    {
        tMaxStat = t;
    }
    
    void setCGTimeAveragingInterval(Mdouble t)
    {
        tIntStat = t;
    }
    
    Mdouble getCGTimeMin()
    {
        return tMinStat;
    }
    
    Mdouble getTimeMaxStat()
    {
        if (std::isnan(tMaxStat))
            return getTimeMax();
        else
            return tMaxStat;
    }
    
    Mdouble getCGTimeAveragingInterval()
    {
        return tIntStat;
    }
    
    bool check_current_time_for_statistics()
    {
        return (getTime() >= getCGTimeMin() && getTime() <= getTimeMaxStat() + getTimeStep());
    }
    
    
    void setCGShape(const char* new_);
    
    void setCGShape(CG_TYPE new_);
    
    CG_TYPE getCGShape()
    {
        return CG_type;
    }
    
    void setN(int nx_, int ny_, int nz_)
    {
        nx = nx_;
        ny = ny_;
        nz = nz_;
    }
    
    void getN(int& nx_, int& ny_, int& nz_)
    {
        nx_ = nx;
        ny_ = ny;
        nz_ = nz;
    }
    
    
    void setCGWidth(Mdouble w)
    {
        setCGWidth2(mathsFunc::square(w));
    }
    
    
    void setCGWidth2(Mdouble new_);
    
    Mdouble getCGWidth()
    {
        return std::sqrt(w2);
    }
    
    Mdouble getCGWidthSquared()
    {
        return w2;
    }
    
    Mdouble getCutoff()
    {
        return cutoff;
    }
    
    Mdouble getCutoff2()
    {
        return mathsFunc::square(cutoff);
    }
    
    
    //~ bool get_boundedDomain() {return boundedDomain;};
    //~ void set_boundedDomain(bool new_) {boundedDomain=new_;};
    
    std::string print_CG();
    
    StatisticsPoint<T> average(std::vector<StatisticsPoint<T> >& P);
    
    virtual void reset_statistics();
    
    void statistics_from_fstat_and_data();
    
    void statistics_from_p3();
    
    void jump_p3c();
    
    void setDoTimeAverage(bool new_)
    {
        doTimeAverage = new_;
    }
    
    bool getDoTimeAverage()
    {
        return doTimeAverage;
    }
    
    
    void setStressTypeForFixedParticles(int new_)
    {
        StressTypeForFixedParticles = new_;
    }
    
    int getStressTypeForFixedParticles()
    {
        return StressTypeForFixedParticles;
    }
    
    
    // to keep compatible with older versions
    void set_infiniteStressForFixedParticles(bool new_)
    {
        StressTypeForFixedParticles = 2 + new_;
    }
    
    
    void setMirrorAtDomainBoundary(Mdouble new_)
    {
        mirrorAtDomainBoundary = new_;
    }
    
    Mdouble getMirrorAtDomainBoundary()
    {
        return mirrorAtDomainBoundary;
    }
    
    
    void setDoVariance(bool new_)
    {
        doVariance = new_;
    }
    
    bool getDoVariance()
    {
        return doVariance;
    }
    
    
    void setDoGradient(bool new_)
    {
        doGradient = new_;
    }
    
    bool getDoGradient()
    {
        return doGradient;
    }
    
    void setSuperExact(bool new_)
    {
        superexact = new_;
    }
    
    bool getSuperExact()
    {
        return superexact;
    }
    
    void setDoIgnoreFixedParticles(bool new_)
    {
        ignoreFixedParticles = new_;
    }
    
    bool getDoIgnoreFixedParticles()
    {
        return ignoreFixedParticles;
    }
    
    void verbose()
    {
        verbosity = 2;
    }
    
    void setVerbosityLevel(int new_)
    {
        verbosity = new_;
    }
    
    int getVerbosityLevel() const
    {
        return verbosity;
    }
    
    void setCGWidthalls(bool new_)
    {
        walls = new_;
    }
    
    bool getCGWidthalls()
    {
        return walls;
    }
    
    void setDoPeriodicWalls(bool new_)
    {
        periodicWalls = new_;
    }
    
    bool getDoPeriodicWalls()
    {
        return periodicWalls;
    }
    
    
    void setCGWidth_over_rmax(Mdouble new_)
    {
        w_over_rmax = new_;
    }
    
    Mdouble getCGWidth_over_rmax()
    {
        return w_over_rmax;
    }
    
    
    void setPositions();
    
    bool readNextDataFile(unsigned int format);
    
    void gather_force_statistics_from_fstat_and_data();
    
    void gather_force_statistics_from_p3c(int version);
    
    void gather_force_statistics_from_p3w(int version, std::vector<int>& index);
    
    void evaluate_force_statistics(int wp = 0);
    
    void evaluate_wall_force_statistics(Vec3D P, int wp = 0);
    
    void jump_fstat();
    
    void initialiseStatistics();
    
    void outputStatistics();
    
    void gatherContactStatistics(unsigned int index1, int index2, Vec3D Contact, Mdouble delta, Mdouble ctheta,
                                 Mdouble fdotn, Mdouble fdott, Vec3D P1_P2_normal_, Vec3D P1_P2_tangential);
    
    void processStatistics(bool usethese);
    
    void finishStatistics();
    
    void write_statistics();
    
    void write_time_average_statistics();
    
    void evaluate_particle_statistics(std::vector<BaseParticle*>::iterator P, int wp = 0);
    
    std::vector<StatisticsPoint<T> > getCGPoints()
    {
        return Points;
    }
    
    Mdouble getXMinStat()
    {
        if (std::isnan(xMinStat))
            return getXMin();
        else
            return xMinStat;
    }
    
    Mdouble getYMinStat()
    {
        if (std::isnan(yMinStat))
            return getYMin();
        else
            return yMinStat;
    }
    
    Mdouble getZMinStat()
    {
        if (std::isnan(zMinStat))
            return getZMin();
        else
            return zMinStat;
    }
    
    Mdouble getXMaxStat()
    {
        if (std::isnan(xMaxStat))
            return getXMax();
        else
            return xMaxStat;
    }
    
    Mdouble getYMaxStat()
    {
        if (std::isnan(yMaxStat))
            return getYMax();
        else
            return yMaxStat;
    }
    
    Mdouble getZMaxStat()
    {
        if (std::isnan(zMaxStat))
            return getZMax();
        else
            return zMaxStat;
    }
    
    void setXMinStat(Mdouble xMinStat_)
    {
        xMinStat = xMinStat_;
    }
    
    void setYMinStat(Mdouble yMinStat_)
    {
        yMinStat = yMinStat_;
    }
    
    void setZMinStat(Mdouble zMinStat_)
    {
        zMinStat = zMinStat_;
    }
    
    void setXMaxStat(Mdouble xMaxStat_)
    {
        xMaxStat = xMaxStat_;
    }
    
    void setYMaxStat(Mdouble yMaxStat_)
    {
        yMaxStat = yMaxStat_;
    }
    
    void setZMaxStat(Mdouble zMaxStat_)
    {
        zMaxStat = zMaxStat_;
    }
    
    int getNTimeAverage()
    {
        return nTimeAverage;
    }
    
    Mdouble setInfinitelyLongDistance();
    
    void set_Polynomial(std::vector<Mdouble> new_coefficients, unsigned int new_dim)
    {
        CGPolynomial.set_polynomial(new_coefficients, new_dim);
    }
    
    void set_Polynomial(Mdouble* new_coefficients, unsigned int num_coeff, unsigned int new_dim)
    {
        CGPolynomial.set_polynomial(new_coefficients, num_coeff, new_dim);
    }
    
    void setPolynomialName(const char* new_name)
    {
        CGPolynomial.setName(new_name);
    }
    
    std::string getPolynomialName()
    {
        return CGPolynomial.getName();
    }
    
    void setDoDoublePoints(bool new_)
    {
        doDoublePoints = new_;
    }
    
    bool getDoDoublePoints()
    {
        return doDoublePoints;
    }
    
    void setNTimeAverageReset(int new_)
    {
        nTimeAverageReset = new_;
    }
    
    bool getNTimeAverageReset()
    {
        return nTimeAverageReset;
    }
    
    void set_rmin(Mdouble new_)
    {
        rmin = new_;
    }
    
    void set_rmax(Mdouble new_)
    {
        rmax = new_;
    }
    
    void set_hmax(Mdouble new_)
    {
        hmax = new_;
    }
    
    Mdouble evaluatePolynomial(Mdouble r)
    {
        return CGPolynomial.evaluate(r);
    }
    
    Mdouble evaluatePolynomialGradient(Mdouble r)
    {
        return CGPolynomial.evaluateGradient(r);
    }
    
    Mdouble evaluateIntegral(Mdouble n1, Mdouble n2, Mdouble t)
    {
        return CGPolynomial.evaluateIntegral(n1, n2, t);
    }
    
    unsigned int getStepSize() const
    {
        return stepSize_;
    }
    
    void setStepSize(unsigned int stepSize)
    {
        stepSize_ = stepSize;
    }
    
    void setIndSpecies(unsigned int indSpecies)
    {
        this->indSpecies = indSpecies;
    }
    //Member Variables
protected:
    
    //General Variables
    StatType statType;
    int nx; 
    int ny; 
    int nz;
    Mdouble xMinStat, xMaxStat, yMinStat, yMaxStat, zMinStat, zMaxStat;
    int nxMirrored, nyMirrored, nzMirrored;
    
    //Storage of points and gradients
    std::vector<StatisticsPoint<T> > Points;
    std::vector<StatisticsPoint<T> > dx;
    std::vector<StatisticsPoint<T> > dy;
    std::vector<StatisticsPoint<T> > dz;
    
    //For time averaging
    std::vector<StatisticsPoint<T> > timeAverage;
    
    std::vector<StatisticsPoint<T> > timeVariance;
    
    std::vector<StatisticsPoint<T> > dxTimeAverage;
    
    std::vector<StatisticsPoint<T> > dyTimeAverage;
    
    std::vector<StatisticsPoint<T> > dzTimeAverage;
    
    bool doTimeAverage;
    
    int nTimeAverageReset;
    
    bool doVariance;
    
    bool doGradient;
    
    int nTimeAverage;
    
    //Coarse graining variables
    CG_TYPE CG_type;
    NORMALIZED_POLYNOMIAL<T> CGPolynomial;
    
    // ///<if true, then the course-graining function will be cut at the domain boundaries and resized to satisfy int(W) = 1
    // bool boundedDomain; 
    
    Mdouble w2;
    
    Mdouble cutoff, cutoff2;
    
    Mdouble w_over_rmax;
    
    //Options that can be set before evaluation
    
    Mdouble tMinStat;
    
    Mdouble tMaxStat;
    
    Mdouble tIntStat;
    
    Mdouble indSpecies;
    
    Mdouble rmin;
    
    Mdouble rmax;
    
    Mdouble hmax;
    
    bool walls;
    
    bool periodicWalls;
    
    bool ignoreFixedParticles;
    
    int StressTypeForFixedParticles;
    
    int verbosity;
    
    int format; 
    
    Mdouble mirrorAtDomainBoundary;
    
    bool isMDCLR;
    bool superexact;
    
    //uses close points to allow calculation of gradients
    bool doDoublePoints;
    
    bool satisfiesInclusionCriteria(BaseParticle* P);
    
    //Variables communicate values between member functions #evaluate_force_statistics and #gatherContactStatistics)used to communicate values between member functions evaluate_force_statistics and gather_force_statistics   
    Vec3D P1;
    Vec3D P2;
    Vec3D P1_P2_normal;
    Mdouble P1_P2_distance;
    Matrix3D P1_P2_NormalStress;
    Matrix3D P1_P2_ContactCoupleStress;
    Vec3D P1_P2_Contact;
    Matrix3D P1_P2_TangentialStress;
    Vec3D P1_P2_NormalTraction;
    Vec3D P1_P2_TangentialTraction;
    MatrixSymmetric3D P1_P2_Fabric;
    Vec3D P1_P2_CollisionalHeatFlux;
    Mdouble P1_P2_Dissipation;
    Mdouble P1_P2_Potential;
    
    bool loadVelocityProfile(const char* filename);
    
    bool loadPositions(const char* filename);
    
    Vec3D getVelocityProfile(Vec3D Position);
    
    std::vector<Vec3D> VelocityProfile;
    
    Vec3D VelocityProfile_Min;
    
    Vec3D VelocityProfile_D;
    
    std::fstream p3p_file;
    std::fstream p3c_file;
    std::fstream p3w_file;
    bool p3Format = true;
    
    bool read_next_from_p3p_file();
    
    void autoSetSystemDimensions();
    
    unsigned int stepSize_;
    
    std::vector<Vec3D> positions_;
    
    bool centerPointInFStat_ = true;
};
 
#include "StatisticsPoint.hcc"
#include "StatisticsVector.hcc"
 
#endif