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