DPMBase.h

Go to the documentation of this file.

//Copyright (c) 2013-2020, 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 DPMBase_H
#define DPMBase_H

//so that the user doesn't have to include string/io manipulations:
#include <string>
#include <iomanip>
// include file for signal handler types
#include <csignal>
//The vector class contains a 3D vector class.
#include "Math/Vector.h"
//This class defines the particle handler
#include "ParticleHandler.h"
//This class defines the base particle (such that not every Driver has to include it)
#include "Particles/BaseParticle.h"
#include "Particles/SphericalParticle.h"
//This class defines the wall handler
#include "WallHandler.h"
//This class defines the boundary handler
#include "BoundaryHandler.h"
#include "PeriodicBoundaryHandler.h"
#include "DomainHandler.h"
//This class defines the interaction handler
#include "InteractionHandler.h"
//This class defines the Species handler
#include "SpeciesHandler.h"
//This class defines the cg handler
#include "CG/CGHandler.h"
//This class defines the possibleContact lists
#ifdef CONTACT_LIST_HGRID
#include "PossibleContactList.h"
#endif
//This class defines the random number generator
#include "Math/RNG.h"
#include "Domain.h"
#include "VTKWriter/ParticleVtkWriter.h"
#include "VTKWriter/WallVTKWriter.h"
#include "VTKWriter/InteractionVTKWriter.h"
#include "VTKWriter/BoundaryVTKWriter.h"
#include "MercuryTime.h"


class DPMBase
{
public:
    
    void constructor();
    
    DPMBase();
    
    DPMBase(const DPMBase& other);
    
    virtual ~DPMBase();
    
    static void incrementRunNumberInFile();
    
    static int readRunNumberFromFile();
    
    void autoNumber();

    std::vector<int> get1DParametersFromRunNumber(int size_x) const;

    std::vector<int> get2DParametersFromRunNumber(int size_x, int size_y) const;

    std::vector<int> get3DParametersFromRunNumber(int size_x, int size_y, int size_z) const;
    
    int launchNewRun(const char* name, bool quick = false);
    
    //setters and getters
    
    void setRunNumber(int runNumber);
    
    int getRunNumber() const;
    
//    void sendParticlesToRoot(unsigned int &numberOfLocalParticles, unsigned int processorId, ParticleHandler &tempParticleHandler) const;
    
    virtual void decompose();
    
    void solve();
    
    virtual void computeOneTimeStep();
    
    void checkSettings();
    
    void forceWriteOutputFiles();
    
    virtual void writeOutputFiles();
    
    void solve(int argc, char* argv[]);
    
    virtual void setupInitialConditions();
    
    virtual void writeXBallsScript() const;
    
    virtual Mdouble getInfo(const BaseParticle& P) const;
    
    ParticleVtkWriter* getVtkWriter() const;
    
    virtual void writeRestartFile();
    
    void writeDataFile();
    
    void writeEneFile();
    
    void writeFStatFile();

    void fillDomainWithParticles(unsigned N=50);

    enum class ReadOptions : int {
        ReadAll,
        ReadNoInteractions,
        ReadNoParticlesAndInteractions
    };

    bool readRestartFile(ReadOptions opt = ReadOptions::ReadAll);

    int readRestartFile(std::string fileName, ReadOptions opt = ReadOptions::ReadAll);

//    /*!
//     * \brief Loads all MD data and plots statistics for all time steps in the .data file
//     */
//    void statisticsFromRestartData(const char *name);
    
    virtual void write(std::ostream& os, bool writeAllParticles = true) const;
    
    virtual void read(std::istream& is, ReadOptions opt = ReadOptions::ReadAll);
    
    virtual BaseWall* readUserDefinedWall(const std::string& type) const
    { return nullptr; }

    virtual void readOld(std::istream& is);

    bool readDataFile(std::string fileName = "", unsigned int format = 0);

    bool readParAndIniFiles(std::string fileName);

    bool readNextDataFile(unsigned int format = 0);

    void readNextFStatFile();

    bool findNextExistingDataFile(Mdouble tMin, bool verbose = true);

    bool readArguments(int argc, char* argv[]);

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

    virtual bool checkParticleForInteraction(const BaseParticle& P);

    virtual bool checkParticleForInteractionLocal(const BaseParticle& P);

    bool checkParticleForInteractionLocalPeriodic(const BaseParticle& P);

    void readSpeciesFromDataFile(bool read = true){readSpeciesFromDataFile_=read;}

    void importParticlesAs(ParticleHandler& particleHandler, InteractionHandler& interactionHandler, const ParticleSpecies* species );

    //getters and setters

    MERCURY_DEPRECATED
    File& getDataFile();

    MERCURY_DEPRECATED
    File& getEneFile();

    MERCURY_DEPRECATED
    File& getFStatFile();

    MERCURY_DEPRECATED
    File& getRestartFile();

    MERCURY_DEPRECATED
    File& getStatFile();

    File& getInteractionFile();

    MERCURY_DEPRECATED
    const File& getDataFile() const;

    MERCURY_DEPRECATED
    const File& getEneFile() const;

    MERCURY_DEPRECATED
    const File& getFStatFile() const;

    MERCURY_DEPRECATED
    const File& getRestartFile() const;

    MERCURY_DEPRECATED
    const File& getStatFile() const;

    const File& getInteractionFile() const;

    const std::string& getName() const;

    void setName(const std::string& name);

    void setName(const char* name);

    void setSaveCount(unsigned int saveCount);

    void setFileType(FileType fileType);

    void setOpenMode(std::fstream::openmode openMode);
    
    
    //other member functions

    void resetFileCounter();

    void closeFiles();

    void setLastSavedTimeStep(unsigned int nextSavedTimeStep);

    Mdouble getTime() const;

    Mdouble getNextTime() const;

    unsigned int getNumberOfTimeSteps() const;

    void setTime(Mdouble time);

    void setTimeMax(Mdouble newTMax);
    
    Mdouble getTimeMax() const;
    
    void setLogarithmicSaveCount(Mdouble logarithmicSaveCountBase);
    
    void setNToWrite(int nToWrite);
    
    int getNToWrite() const;

#ifdef CONTACT_LIST_HGRID

        PossibleContactList& getPossibleContactList();
#endif
    
    void setRotation(bool rotation)
    { rotation_ = rotation; }

    bool getRotation() const
    { return rotation_; }

    void setWallsWriteVTK(FileType writeWallsVTK);

    void setWallsWriteVTK(bool);

    void setInteractionsWriteVTK(bool);

    void setParticlesWriteVTK(bool writeParticlesVTK);

    void setSuperquadricParticlesWriteVTK(bool writeSuperquadricParticlesVTK);

    FileType getWallsWriteVTK() const;

    bool getParticlesWriteVTK() const;

    bool getSuperquadricParticlesWriteVTK() const;

    Mdouble getXMin() const
    { return min_.x(); }

    Mdouble getXMax() const
    { return max_.x(); }

    Mdouble getYMin() const
    { return min_.y(); }

    Mdouble getYMax() const
    { return max_.y(); }

    Mdouble getZMin() const
    { return min_.z(); }

    Mdouble getZMax() const
    { return max_.z(); }

    /*
     * \brief Returns the minimum coordinates of the problem domain.
     */
    Vec3D getMin() const
    { return min_; }

    /*
     * \brief Returns the maximum coordinates of the problem domain.
     */
    Vec3D getMax() const
    { return max_; }

    void setXMin(Mdouble newXMin);

    void setYMin(Mdouble newYMin);

    void setZMin(Mdouble newZMin);

    void setXMax(Mdouble newXMax);

    void setYMax(Mdouble newYMax);

    void setZMax(Mdouble newZMax);

    void setMax(const Vec3D& max);

    void setMax(Mdouble, Mdouble, Mdouble);

    void setDomain(const Vec3D& min, const Vec3D& max);

    void setMin(const Vec3D& min);

    void setMin(Mdouble, Mdouble, Mdouble);


    void setTimeStep(Mdouble newDt);

    Mdouble getTimeStep() const;

    /* Sets the number of omp threads */
    void setNumberOfOMPThreads(int numberOfOMPThreads);

    /* Returns the number of omp threads */
    int getNumberOfOMPThreads() const;

    void setXBallsColourMode(int newCMode);

    int getXBallsColourMode() const;

    void setXBallsVectorScale(double newVScale);

    double getXBallsVectorScale() const;

    void setXBallsAdditionalArguments(std::string newXBArgs);

    std::string getXBallsAdditionalArguments() const;

    void setXBallsScale(Mdouble newScale);

    double getXBallsScale() const;

    void setGravity(Vec3D newGravity);

    Vec3D getGravity() const;

    void setDimension(unsigned int newDim);

    void setSystemDimensions(unsigned int newDim);

    unsigned int getSystemDimensions() const;

    void setParticleDimensions(unsigned int particleDimensions);

    unsigned int getParticleDimensions() const;

    std::string getRestartVersion() const;

    void setRestartVersion(std::string newRV);

    bool getRestarted() const;

    void setRestarted(bool newRestartedFlag);

    bool getAppend() const;

    void setAppend(bool newAppendFlag);

    Mdouble getElasticEnergy() const;

    Mdouble getKineticEnergy() const;

    Mdouble getGravitationalEnergy() const;

    Mdouble getRotationalEnergy() const;

    Mdouble getTotalEnergy() const;

    Mdouble getTotalMass() const;

    Vec3D getCentreOfMass() const;

    Vec3D getTotalMomentum() const;

    double getCPUTime() { return clock_.getCPUTime(); }
    
    double getWallTime() { return clock_.getWallTime(); }
    
    
    static bool areInContact(const BaseParticle* pI, const BaseParticle* pJ);


    virtual void hGridInsertParticle(BaseParticle* obj UNUSED);

    virtual void hGridUpdateParticle(BaseParticle* obj UNUSED);

    virtual void hGridRemoveParticle(BaseParticle* obj UNUSED);

    virtual void hGridUpdateMove(BaseParticle*, Mdouble);

    bool mpiIsInCommunicationZone(BaseParticle* particle);

    bool mpiInsertParticleCheck(BaseParticle* P);

    void insertGhostParticle(BaseParticle* P);

    void updateGhostGrid(BaseParticle* P);

    virtual 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 setNumberOfDomains(std::vector<unsigned> direction);

    enum class DomainSplit {X, Y, Z, XY, XZ, YZ, XYZ};

    void splitDomain(DomainSplit domainSplit);

    std::vector<unsigned> getNumberOfDomains();

    Domain* getCurrentDomain();

    void removeOldFiles() const;

    virtual void hGridGetInteractingParticleList(BaseParticle* obj, std::vector<BaseParticle*>& list)
    {};

    virtual void computeWallForces(BaseWall* w);

    virtual bool getHGridUpdateEachTimeStep() const;

    void setMeanVelocity(Vec3D V_mean_goal);

    void setMeanVelocityAndKineticEnergy(Vec3D V_mean_goal, Mdouble Ek_goal);

    Mdouble getTotalVolume() const;

    Matrix3D getKineticStress() const;

    Matrix3D getStaticStress() const;

    Matrix3D getTotalStress() const;
    
    virtual void handleParticleRemoval(unsigned int id);
    
    virtual void handleParticleAddition(unsigned int id, BaseParticle* p);

    //functions that should only be used in the class definitions
protected:

    virtual void computeAllForces();

    virtual void computeInternalForces(BaseParticle*);

    virtual void computeInternalForce(BaseParticle*, BaseParticle*);

    virtual void computeExternalForces(BaseParticle*);

    void computeForcesDueToWalls(BaseParticle*, BaseWall*);

    virtual void actionsOnRestart();

    virtual void actionsBeforeTimeLoop();

    virtual void hGridActionsBeforeTimeLoop();

    virtual void hGridActionsBeforeTimeStep();

    virtual void actionsBeforeTimeStep();

    virtual void computeAdditionalForces() {}
    
    virtual void actionsAfterSolve();

    virtual void actionsAfterTimeStep();

    void writeVTKFiles() const;

    virtual void outputXBallsData(std::ostream& os) const;

    virtual void outputXBallsDataParticle(unsigned int i, unsigned int format, std::ostream& os) const;

    virtual void writeEneHeader(std::ostream& os) const;

    virtual void writeFstatHeader(std::ostream& os) const;

    virtual void writeEneTimeStep(std::ostream& os) const;

    // Functions for statistics
    virtual void initialiseStatistics();

    virtual void outputStatistics();

    void gatherContactStatistics();

    virtual void processStatistics(bool);

    virtual void finishStatistics();

    virtual void integrateBeforeForceComputation();

    virtual void integrateAfterForceComputation();

    virtual void checkInteractionWithBoundaries();

    virtual void hGridActionsBeforeIntegration();

    virtual void hGridActionsAfterIntegration();

    void setFixedParticles(unsigned int n);

    virtual void printTime() const;

    virtual bool continueSolve() const;

    void outputInteractionDetails() const;

    bool isTimeEqualTo(Mdouble time) const;

    void removeDuplicatePeriodicParticles();

    void checkAndDuplicatePeriodicParticles();

    void performGhostParticleUpdate();

    void deleteGhostParticles(std::set<BaseParticle*>& particlesToBeDeleted);

    void synchroniseParticle(BaseParticle*, unsigned fromProcessor = 0);

    void performGhostVelocityUpdate();


    static void signalHandler(int signal);

    void setSoftStop();

private:

    static volatile sig_atomic_t continueFlag_;

    /*The number of openmp (symmetric multiprocessing threads)*/
    int numberOfOMPThreads_;

    unsigned int systemDimensions_;

    unsigned int particleDimensions_;

    Vec3D gravity_;

    std::vector<unsigned> numberOfDomains_;

    Vec3D min_;
    Vec3D max_;

    Mdouble time_;

    unsigned int numberOfTimeSteps_;

    Mdouble timeStep_;

    Mdouble timeMax_;

    std::string restartVersion_;

    bool restarted_;

    bool append_;

    bool rotation_;

    FileType writeWallsVTK_;

    bool writeParticlesVTK_;

    bool writeSuperquadricParticlesVTK_;

    ParticleVtkWriter* vtkWriter_;

    WallVTKWriter wallVTKWriter_;

    InteractionVTKWriter interactionVTKWriter_;

    BoundaryVTKWriter boundaryVTKWriter_;

    //This is the private data that is only used by the xballs output

    int xBallsColourMode_;

    Mdouble xBallsVectorScale_;

    Mdouble xBallsScale_;

    std::string xBallsAdditionalArguments_;

    int runNumber_;

    std::string name_;
    
    // defines a Macro for creating an instance of class PossibleContactList. See PossibleContactList.h
#ifdef CONTACT_LIST_HGRID
    PossibleContactList possibleContactList;
#endif
    
    bool readSpeciesFromDataFile_;
    
    int nToWrite_;

public:
    SpeciesHandler speciesHandler;
    
    RNG random;

    ParticleHandler particleHandler;

    ParticleHandler paoloParticleHandler;

    WallHandler wallHandler;

    BoundaryHandler boundaryHandler;

    PeriodicBoundaryHandler periodicBoundaryHandler;

    DomainHandler domainHandler;

    InteractionHandler interactionHandler;


    CGHandler cgHandler;

    File dataFile;

    File fStatFile;

    File eneFile;

    File restartFile;

    File statFile;

    File interactionFile;

    Time clock_;
    
    void writePythonFileForVTKVisualisation() const;
};

#endif