StatisticsVector.h

Go to the documentation of this file.

//Copyright (c) 2013-2014, The MercuryDPM Developers Team. All rights reserved.
//For the list of developers, see <http://www.MercuryDPM.org/Team>.
//
//Redistribution and use in source and binary forms, with or without
//modification, are permitted provided that the following conditions are met:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name MercuryDPM nor the
//    names of its contributors may be used to endorse or promote products
//    derived from this software without specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
//ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
//WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
//DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
//LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
//ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
//(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
//SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#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)
            std::cerr << "Warning in set_nz: dimension less than 3" << std::endl;
    }
    
    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 new_);

    CG 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 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 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 read_next_from_p3p_file();

    void auto_setdim();

    unsigned int stepSize_;

    std::vector<Vec3D> positions_;
};

#include "StatisticsPoint.hcc"
#include "StatisticsVector.hcc"

#endif