DPMBase.cc

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.


#include "DPMBase.h"
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <cmath>
#include <fstream>
#include <cstdlib>
#include <limits>
#include <string>
#include <sstream>
#include <cstdio>
#include <cstring>
//This is only used to change the file permission of the xball script create, at some point this code may be moved from this file to a different file.
#include <sys/types.h>
#include <sys/stat.h>
#include <Species/LinearViscoelasticSlidingFrictionSpecies.h>
#include "Interactions/Interaction.h"

#include "Species/Species.h"
#include "Species/LinearViscoelasticSpecies.h"
#include "Species/FrictionForceSpecies/SlidingFrictionSpecies.h"
#include "Particles/BaseParticle.h"

//This is part of this class and just separates out the stuff to do with xballs.
#include "CMakeDefinitions.h"
#include "DPMBaseXBalls.icc"
#include "Logger.h"
#include "Particles/BaseParticle.h"
#include "Walls/BaseWall.h"
#include "Walls/InfiniteWall.h"
#include "Boundaries/PeriodicBoundary.h"

void logWriteAndDie(std::string module, std::string message)
{
    std::cerr << "A fatal   error has occured"
              << "\n Module  :" << module
              << "\n Message :" << message << std::endl;
    
    std::exit(-1);
}

std::ostream& operator<<(std::ostream& os, DPMBase& md)
{
    md.write(os);
    return os;
}

DPMBase::DPMBase(const FilesAndRunNumber& other)
        : FilesAndRunNumber(other)
{
    constructor();
}

DPMBase::DPMBase(const DPMBase& other)
{
    elasticEnergy_ = other.elasticEnergy_;
    systemDimensions_ = other.systemDimensions_;
    particleDimensions_ = other.particleDimensions_;
    gravity_ = other.gravity_;
    xMin_ = other.xMin_;
    xMax_ = other.xMax_;
    yMin_ = other.yMin_;
    yMax_ = other.yMax_;
    zMin_ = other.zMin_;
    zMax_ = other.zMax_;
    time_ = other.time_;
    timeStep_ = other.timeStep_;
    ntimeSteps_ = other.ntimeSteps_;
    timeMax_ = other.timeMax_;
    restartVersion_ = other.restartVersion_; //to read new and old restart data
    restarted_ = other.restarted_; //to see if it was restarted or not
    append_ = other.append_;
    rotation_ = other.rotation_;
    xBallsColourMode_ = other.xBallsColourMode_; // sets the xballs argument cmode (see xballs.txt)
    xBallsVectorScale_ = other.xBallsVectorScale_; // sets the xballs argument vscale (see xballs.txt)
    xBallsScale_ = other.xBallsScale_; // sets the xballs argument scale (see xballs.txt)
    xBallsAdditionalArguments_ = other.xBallsAdditionalArguments_; // std::string where additional xballs argument can be specified (see xballs.txt)
    writeWallsVTK_ = other.writeWallsVTK_;
    writeParticlesVTK_ = other.writeParticlesVTK_;
#ifdef CONTACT_LIST_HGRID
    possibleContactList=other.possibleContactList;
#endif
    random = other.random;

    boundaryHandler.setDPMBase(this);
    particleHandler.setDPMBase(this);
    interactionHandler.setDPMBase(this);
    speciesHandler.setDPMBase(this);
    wallHandler.setDPMBase(this);

    speciesHandler = other.speciesHandler;
    particleHandler = other.particleHandler;
    wallHandler = other.wallHandler;
    boundaryHandler = other.boundaryHandler;
    interactionHandler = other.interactionHandler;
}

DPMBase::DPMBase()
{
    constructor();
}

DPMBase::~DPMBase()
{

}

void DPMBase::solve(int argc, char* argv[])
{
    readArguments(argc, argv);
    solve();
}

Mdouble DPMBase::getTime() const
{
    return time_;
}

unsigned int DPMBase::getNtimeSteps() const
{
    return ntimeSteps_;
}

void DPMBase::setTime(Mdouble time)
{
    Mdouble diff = time_ - time;
    time_ = time;
    //this sets the interaction timestamp, so each interaction has the right time
    for (auto i : interactionHandler)
    {
        i->setTimeStamp(i->getTimeStamp() - diff);
    }
}

void DPMBase::setTimeMax(Mdouble newTMax)
{
    if (newTMax >= 0)
    {
        timeMax_ = newTMax;
    }
    else
    {
        logger(ERROR, "Error in setTimeMax, newTMax=%", newTMax);
    }
}

Mdouble DPMBase::getTimeMax() const
{
    return timeMax_;
}
#ifdef CONTACT_LIST_HGRID
PossibleContactList& DPMBase::getPossibleContactList()
{   
    return possibleContactList;
}
#endif

void DPMBase::setRotation(bool newRotFlag)
{
    rotation_ = newRotFlag;
}

bool DPMBase::getRotation() const
{
    return rotation_;
}

void DPMBase::setWallsWriteVTK(FileType writeWallsVTK)
{
    writeWallsVTK_ = writeWallsVTK;
}

void DPMBase::setParticlesWriteVTK(bool writeParticlesVTK)
{
    writeParticlesVTK_ = writeParticlesVTK;
}

FileType DPMBase::getWallsWriteVTK() const
{
    return writeWallsVTK_;
}

bool DPMBase::getParticlesWriteVTK() const
{
    return writeParticlesVTK_;
}

Mdouble DPMBase::getXMin() const
{
    return xMin_;
}

Mdouble DPMBase::getXMax() const
{
    return xMax_;
}

Mdouble DPMBase::getYMin() const
{
    return yMin_;
}

Mdouble DPMBase::getYMax() const
{
    return yMax_;
}

Mdouble DPMBase::getZMin() const
{
    return zMin_;
}

Mdouble DPMBase::getZMax() const
{
    return zMax_;
}

int DPMBase::getAddLarge() const
{
    return addL_;
}

int DPMBase::getAddSmall() const
{
    return addS_;
}

void DPMBase::setAddSmall(int addS)
{
    addS_ = addS;
}

void DPMBase::setAddLarge(int addL)
{
    addL_ = addL;
}

void DPMBase::setXMin(Mdouble xMin)
{
    if (xMin > getXMax())
    {
        logger(WARN, "Warning in setXMin(%): value cannot be larger than xMax=%", xMin, getXMax());
    }
    xMin_ = xMin;
}

void DPMBase::setYMin(Mdouble yMin)
{
    if (yMin > getYMax())
    {
        logger(WARN, "Warning in setYMin(%): value cannot be larger than yMax=%", yMin, getYMax());
    }
    yMin_ = yMin;
}

void DPMBase::setZMin(Mdouble zMin)
{
    if (zMin > getZMax())
    {
        logger(WARN, "Warning in setZMin(%): value cannot be larger than zMax=%", zMin, getZMax());
    }
    zMin_ = zMin;
}

void DPMBase::setMax(Vec3D max)
{
    setXMax(max.X);
    setYMax(max.Y);
    setZMax(max.Z);
}

void DPMBase::setMin(Vec3D min)
{
    setXMin(min.X);
    setYMin(min.Y);
    setZMin(min.Z);
}

void DPMBase::setXMax(Mdouble xMax)
{
    if (xMax < getXMin())
    {
        logger(WARN, "Warning in setXMax(%): value cannot be smaller than xMin=%", xMax, getXMin());
    }
    xMax_ = xMax;
}

void DPMBase::setYMax(Mdouble yMax)
{
    if (yMax < getYMin())
    {
        logger(WARN, "Warning in setYMax(%): value cannot be smaller than yMin=%", yMax, getYMin());
    }
    yMax_ = yMax;
}

void DPMBase::setZMax(Mdouble zMax)
{
    if (zMax < getZMin())
    {
        logger(WARN, "Warning in setZMax(%): value cannot be smaller than zMin=%", zMax, getZMin());
    }
    zMax_ = zMax;
}

void DPMBase::setTimeStep(Mdouble timeStep)
{
    if (timeStep >= 0.0)
    {
        timeStep_ = timeStep;
    } else
    {
        logger(ERROR, "Error in setTimeStep, proposed time step %", timeStep);
    }
}

Mdouble DPMBase::getTimeStep() const
{
    return timeStep_;
}

void DPMBase::setXBallsColourMode(int newCMode)
{
    xBallsColourMode_ = newCMode;
}

int DPMBase::getXBallsColourMode() const
{
    return xBallsColourMode_;
}

void DPMBase::setXBallsVectorScale(double newVScale)
{
    xBallsVectorScale_ = newVScale;
}

double DPMBase::getXBallsVectorScale() const
{
    return xBallsVectorScale_;
}

void DPMBase::setXBallsAdditionalArguments(std::string newXBArgs)
{
    xBallsAdditionalArguments_ = newXBArgs.c_str();
}

std::string DPMBase::getXBallsAdditionalArguments() const
{
    return xBallsAdditionalArguments_;
}

void DPMBase::setXBallsScale(Mdouble newScale)
{
    xBallsScale_ = newScale;
}

double DPMBase::getXBallsScale() const
{
    return xBallsScale_;
}

void DPMBase::setGravity(Vec3D newGravity)
{
    gravity_ = newGravity;
}

Vec3D DPMBase::getGravity() const
{
    return gravity_;
}

void DPMBase::setDimension(unsigned int newDim)
{
    setSystemDimensions(newDim);
    setParticleDimensions(newDim);
}

void DPMBase::setSystemDimensions(unsigned int newDim)
{
    if (newDim >= 1 && newDim <= 3)
        systemDimensions_ = newDim;
    else
    {
        logger(ERROR,"Error in setSystemDimensions, input argument is %", newDim);
    }
}

unsigned int DPMBase::getSystemDimensions() const
{
    return systemDimensions_;
}

void DPMBase::setParticleDimensions(unsigned int particleDimensions)
{
    if (particleDimensions >= 1 && particleDimensions <= 3)
    {
        particleDimensions_ = particleDimensions;
        particleHandler.computeAllMasses();
    } else
    {
        logger(WARN, "Error in setParticleDimensions, input argument is %", particleDimensions);
    }
}

unsigned int DPMBase::getParticleDimensions() const
{
    return particleDimensions_;
}

std::string DPMBase::getRestartVersion() const
{
    return restartVersion_;
}

void DPMBase::setRestartVersion(std::string newRV)
{
    restartVersion_ = newRV;
}

bool DPMBase::getRestarted() const
{
    return restarted_;
}

void DPMBase::setRestarted(bool newRestartedFlag)
{
    restarted_ = newRestartedFlag;
    //setAppend(new_);
}

bool DPMBase::getAppend() const
{
    return append_;
}

void DPMBase::setAppend(bool newAppendFlag)
{
    append_ = newAppendFlag;
}

Mdouble DPMBase::getElasticEnergy() const
{
    Mdouble elasticEnergy = 0.0;
    for (const BaseInteraction* c : interactionHandler)
    {
        elasticEnergy += c->getElasticEnergy();
    }
    return elasticEnergy;
}

Mdouble DPMBase::getKineticEnergy() const
{
    Mdouble kineticEnergy = 0;
    for (const BaseParticle* const p : particleHandler)
    {
        if (!(p->isFixed()))
        {
            kineticEnergy += .5 * p->getMass() * p->getVelocity().getLengthSquared();
        }
    }
    return kineticEnergy;
}

//finding the global **potential** energy
Mdouble DPMBase::getGravitationalEnergy() const
{
    Mdouble gravitationalEnergy = 0;
    for (const BaseParticle* const p : particleHandler)
    {
        if (!(p->isFixed()))
        {
            gravitationalEnergy += p->getMass() * Vec3D::dot((-getGravity()), p->getPosition());
        }
    }
    return gravitationalEnergy;

}

double DPMBase::getInfo(const BaseParticle& p) const
{
    return p.getSpecies()->getId(); // was getIndex()
}

bool DPMBase::areInContact(const BaseParticle* pI, const BaseParticle* pJ) const
{
    return (pI != pJ && Vec3D::getDistanceSquared(pI->getPosition(), pJ->getPosition()) <
                       mathsFunc::square(pI->getInteractionRadius() + pJ->getInteractionRadius()));
}

void DPMBase::actionsBeforeTimeLoop()
{
}

void DPMBase::hGridActionsBeforeTimeLoop()
{
}

void DPMBase::actionsOnRestart()
{
}

void DPMBase::hGridActionsBeforeTimeStep()
{
}

void DPMBase::hGridInsertParticle(BaseParticle* obj UNUSED)
{
}

void DPMBase::hGridUpdateParticle(BaseParticle* obj UNUSED)
{
}

void DPMBase::hGridRemoveParticle(BaseParticle* obj UNUSED)
{
}

bool DPMBase::getHGridUpdateEachTimeStep() const
{
    return true;
}

void DPMBase::actionsBeforeTimeStep()
{
}

void DPMBase::actionsAfterSolve()
{
}

void DPMBase::actionsAfterTimeStep()
{
}

void DPMBase::initialiseStatistics()
{
}

void DPMBase::outputStatistics()
{
}

void DPMBase::gatherContactStatistics()
{
    for (BaseInteraction* c : interactionHandler)
    {
        c->gatherContactStatistics();
    }
}

void DPMBase::gatherContactStatistics(unsigned int index1 UNUSED, int index2 UNUSED, Vec3D Contact UNUSED,
                                      Mdouble delta UNUSED, Mdouble ctheta UNUSED, Mdouble fdotn UNUSED,
                                      Mdouble fdott UNUSED, Vec3D P1_P2_normal_ UNUSED, Vec3D P1_P2_tangential UNUSED)
{
}

void DPMBase::processStatistics(bool usethese UNUSED)
{
}

void DPMBase::finishStatistics()
{
}

void DPMBase::hGridUpdateMove(BaseParticle*, Mdouble)
{
}

void DPMBase::hGridActionsBeforeIntegration()
{
}

void DPMBase::hGridActionsAfterIntegration()
{
}

void DPMBase::broadPhase(BaseParticle* i)
{
    for (auto it = particleHandler.begin(); (*it) != i; ++it)
    {
        computeInternalForces(i, *it);
    }
}

void DPMBase::setFixedParticles(unsigned int n)
{
    for (unsigned int i = 0; i < std::min(particleHandler.getNumberOfObjects(), n); i++)
        particleHandler.getObject(i)->fixParticle();
}

void DPMBase::printTime() const
{
    std::cout << "t=" << std::setprecision(3) << std::left << std::setw(6) << getTime()
              << ", tmax=" << std::setprecision(3) << std::left << std::setw(6) << getTimeMax()
              << std::endl;
    std::cout.flush();
}

bool DPMBase::continueSolve() const
{
    return true;
}

void DPMBase::constructor()
{
    //These are the particle parameters like dissipation etc..
    setRestarted(false);

    //This sets the maximum number of particles
    boundaryHandler.setDPMBase(this);
    speciesHandler.setDPMBase(this);
    particleHandler.setDPMBase(this);
    interactionHandler.setDPMBase(this);
    wallHandler.setDPMBase(this);
    interactionHandler.setDPMBase(this);

    particleHandler.setStorageCapacity(2);
    setSystemDimensions(2);
    setParticleDimensions(2);

    // Set gravitational acceleration
    gravity_ = Vec3D(0.0, -9.8, 0.0);

    //This is the parameter of the numerical part
    time_ = 0.0;
    timeStep_ = 0.0; // if dt is not user-specified, this is set in actionsBeforeTimeLoop()
    ntimeSteps_ = 0;
    setSaveCount(20);
    timeMax_ = 1.0;

    //This sets the default xballs domain
    xMin_ = 0.0;
    xMax_ = 0.01;
    yMin_ = 0.0;
    yMax_ = 0.01;
    zMin_ = 0.0;
    zMax_ = 0.0;

    //sets the default write particles data in VTK format flag to false
    writeParticlesVTK_ = false;
    writeWallsVTK_ = FileType::NO_FILE;

    loggerOutput->onFatal = logWriteAndDie;

    setName("");

    //Default mode is energy with no scale of the vectors
    xBallsColourMode_ = 0;
    xBallsVectorScale_ = -1;
    xBallsScale_ = -1;
    xBallsAdditionalArguments_ = "";
    setAppend(false);
//TODO
    //The default random seed is 0
    random.setRandomSeed(0);

#ifdef DEBUG_OUTPUT
    std::cerr << "MD problem constructor finished " << std::endl;
#endif
}

void DPMBase::setupInitialConditions()
{
}

void DPMBase::writeEneHeader(std::ostream& os) const
{
    //only write if we don't restart
    if (getAppend())
        return;
    
    long width = os.precision() + 6;
    os << std::setw(width) << "t" << " " << std::setw(width)
       << "ene_gra" << " " << std::setw(width)
       << "ene_kin" << " " << std::setw(width)
       << "ene_rot" << " " << std::setw(width)
       << "ene_ela" << " " << std::setw(width)
       << "X_COM" << " " << std::setw(width)
       << "Y_COM" << " " << std::setw(width)
       << "Z_COM" << std::endl;
}

void DPMBase::writeFstatHeader(std::ostream& os) const
{

    // line #1: time, volume fraction
    // line #2: wall box: wx0, wy0, wz0, wx1, wy1, wz1
    // line #3: radii-min-max & moments: rad_min, rad_max, r1, r2, r3, r4
    os << "#"
       << " " << getTime()
       << " " << 0
       << std::endl;
    os << "#"
       << " " << getXMin()
       << " " << getYMin()
       << " " << getZMin()
       << " " << getXMax()
       << " " << getYMax()
       << " " << getZMax()
       << std::endl;
    os << "#"
       << " ";
    if (particleHandler.getSmallestParticle())
    {
        os << particleHandler.getSmallestParticle()->getRadius();
    } else
    {
        os << std::numeric_limits<double>::quiet_NaN();
    }
    os << " ";
    if (particleHandler.getLargestParticle())
    {
        os << particleHandler.getLargestParticle()->getRadius();
    } else
    {
        os << std::numeric_limits<double>::quiet_NaN();
    }
    os << " " << 0
       << " " << 0
       << " " << 0
       << " " << 0
       << std::endl;
    //B: write data
    for (BaseInteraction* c : interactionHandler)
    {
        c->writeToFStat(os);
    }
}

void DPMBase::writeEneTimestep(std::ostream& os) const
{
    Mdouble ene_kin = 0, ene_elastic = 0, ene_rot = 0, ene_gra = 0, mass_sum = 0, x_masslength = 0, y_masslength = 0, z_masslength = 0;

    for (BaseParticle* p : particleHandler)
        if (!p->isFixed())
        {
            ene_kin += .5 * p->getMass() * p->getVelocity().getLengthSquared();
            ene_rot += .5 * p->getInertia() * p->getAngularVelocity().getLengthSquared();
            ene_gra -= p->getMass() * Vec3D::dot(getGravity(), p->getPosition());
            mass_sum += p->getMass();
            x_masslength += p->getMass() * p->getPosition().X;
            y_masslength += p->getMass() * p->getPosition().Y;
            z_masslength += p->getMass() * p->getPosition().Z;
        } //end for loop over Particles

    ene_elastic = getElasticEnergy();

    long width = os.precision() + 6;
    os << std::setw(width) << getTime()
       << " " << std::setw(width) << ene_gra
       << " " << std::setw(width) << ene_kin
       << " " << std::setw(width) << ene_rot
       << " " << std::setw(width) << ene_elastic
       << " " << std::setw(width)
       << (mass_sum != 0.0 ? x_masslength / mass_sum : std::numeric_limits<double>::quiet_NaN())
       << " " << std::setw(width)
       << (mass_sum != 0.0 ? y_masslength / mass_sum : std::numeric_limits<double>::quiet_NaN())
       << " " << std::setw(width)
       << (mass_sum != 0.0 ? z_masslength / mass_sum : std::numeric_limits<double>::quiet_NaN()) << std::endl;

    //sliding = sticking = 0;
}


void DPMBase::writeVTK() const {
    if ((getWallsWriteVTK() == FileType::ONE_FILE && getTime() == 0.0)
        || getWallsWriteVTK() == FileType::MULTIPLE_FILES
        || getWallsWriteVTK() == FileType::MULTIPLE_FILES_PADDED)
    {
        wallHandler.writeVTK();
    }
    if (getParticlesWriteVTK() == true)
    {
        particleHandler.writeVTK();
    }
    interactionHandler.writeVTK();

    if (getTime()==0 && interactionHandler.getWriteVTK() != FileType::NO_FILE)
    {
        logger(INFO,"Writing python script for paraview visualisation");
        helpers::writeToFile(getName()+".paraview.py",
                             "#script to visualise the output of data2pvd of MercuryDPM in paraview.\n"
                              "#usage: change the path below to your own path, open paraview\n"
                              "#Tools->Python Shell->Run Script->VisualisationScript.py\n"
                              "#or run paraview --script=VisualisationScript.py \n"
                              "\n"
                              "try: paraview.simple\n"
                              "except: from paraview.simple import *\n"
                              "import glob\n"
                              //"import os\n"
                              //"os.chdir(\"/Users/weinhartt/Mercury/Alpha/cmake-build-debug/Drivers/USER/Thomas/Mainz\")\n"
                              "\n"
                              "particles = XMLUnstructuredGridReader(FileName=glob.glob('./"
                             + getName() + "Particle_*.vtu'))\n"
                              "\n"
                              "walls = XMLUnstructuredGridReader(FileName=glob.glob('./"
                             + getName() + "Wall_*.vtu'))\n"
                              "\n"
                              "interactions = XMLUnstructuredGridReader(FileName=glob.glob('./"
                             + getName() + "Interaction_*.vtu'))\n"
                              "\n"
                              "#Make the glyph, see visualisation guide on the docs\n"
                              "glyphP = Glyph(particles)\n"
                              "glyphP.GlyphType = 'Sphere'\n"
                              "glyphP.Scalars = 'Radius'\n"
                              "glyphP.Vectors = 'None'\n"
                              "glyphP.ScaleMode = 'scalar'\n"
                              "glyphP.ScaleFactor = 2\n"
                              "glyphP.GlyphMode = 'All Points'\n"
                              "\n"
                              "#Make the glyph, see visualisation guide on the docs\n"
                              "glyphI = Glyph(interactions)\n"
                              "glyphI.GlyphType = 'Sphere'\n"
                              "glyphI.Scalars = 'ContactRadius'\n"
                              "glyphI.Vectors = 'None'\n"
                              "glyphI.ScaleMode = 'scalar'\n"
                              "glyphI.ScaleFactor = 10\n" //5 times too large
                              "glyphI.GlyphMode = 'All Points'\n"
                              "\n"
                              "#show on screen\n"
                              "Show(glyphI)\n"
                              "Show(glyphP)\n"
                              "Show(walls)\n"
                              "Render()\n"
                              "ResetCamera()");
    }
}

void DPMBase::outputXBallsData(std::ostream& os) const
{
    //Set the correct formation based of dimension if the formation is not specified by the user
    unsigned int format;
    switch (getSystemDimensions())
    {
        case 2:
            format = 8;
            break;
        case 3:
            format = 14;
            break;
        default:
            std::cerr << "Unknown systemdimension" << std::endl;
            exit(-1);
    }

    // This outputs the location of walls and how many particles there are to file this is required by the xballs plotting
    if (format != 14) // dim = 1 or 2
    {
        os << particleHandler.getNumberOfObjects()
           << " " << getTime()
           << " " << getXMin()
           << " " << getYMin()
           << " " << getXMax()
           << " " << getYMax()
           << " " << std::endl;
    } else
    {
        //dim==3
        os << particleHandler.getNumberOfObjects()
           << " " << getTime()
           << " " << getXMin()
           << " " << getYMin()
           << " " << getZMin()
           << " " << getXMax()
           << " " << getYMax()
           << " " << getZMax()
           << " " << std::endl;
    }
    // This outputs the particle data
    for (unsigned int i = 0; i < particleHandler.getNumberOfObjects(); i++)
    {
        outputXBallsDataParticle(i, format, os);
    }
#ifdef DEBUG_OUTPUT
    std::cerr << "Have output the properties of the problem to disk " << std::endl;
#endif
}

bool DPMBase::readDataFile(const std::string fileName, unsigned int format)
{
    std::string oldFileName = dataFile.getName();
    //This opens the file the data will be recalled from
    dataFile.setName(fileName);
    dataFile.open(std::fstream::in);
    if (!dataFile.getFstream().is_open() || dataFile.getFstream().bad())
    {
        std::cout << "Loading data file " << fileName << " failed" << std::endl;
        return false;
    }

    //dataFile.getFileType() is ignored
    FileType fileTypeData = dataFile.getFileType();
    dataFile.setFileType(FileType::ONE_FILE);
    readNextDataFile(format);
    dataFile.setFileType(fileTypeData);

    dataFile.close();
    dataFile.setName(oldFileName);

    //std::cout << "Loaded data file " << filename << " (t=" << getTime() << ")" << std::endl;
    return true;
}

bool DPMBase::readParAndIniFiles(const std::string fileName)
{
    //Opens the par.ini file
    std::fstream file;
    file.open(fileName, std::fstream::in);
    if (!file.is_open() || file.bad())
    {
        //std::cout << "Loading par.ini file " << filename << " failed" << std::endl;
        return false;
    }
    
    Mdouble doubleValue;
    int integerValue;
    
    // inputfile par.ini 
    // line 1 =============================================================
    // Example: 1 1 0
    //   1: integer (0|1) switches from non-periodic to periodic
    //      integer (5|6) does 2D integration only (y-coordinates fixed)
    //                    and switches from non-periodic to periodic
    //      integer (11) uses a quarter system with circular b.c.
    file >> integerValue;
    //~ std::cout << "11" << integerValue << std::endl;
    if (integerValue == 0)
    {
        InfiniteWall w0;
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(-1, 0, 0), Vec3D(getXMin(), 0, 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(1, 0, 0), Vec3D(getXMax(), 0, 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0, -1, 0), Vec3D(0, getYMin(), 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0, 1, 0), Vec3D(0, getYMax(), 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0, 0, -1), Vec3D(0, 0, getZMin()));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0, 0, 1), Vec3D(0, 0, getZMax()));
        wallHandler.copyAndAddObject(w0);
    } else if (integerValue == 1)
    {
        PeriodicBoundary b0;
        b0.set(Vec3D(1, 0, 0), getXMin(), getXMax());
        boundaryHandler.copyAndAddObject(b0);
        b0.set(Vec3D(0, 1, 0), getYMin(), getYMax());
        boundaryHandler.copyAndAddObject(b0);
        b0.set(Vec3D(0, 0, 1), getZMin(), getZMax());
        boundaryHandler.copyAndAddObject(b0);
    } else if (integerValue == 5)
    {
        InfiniteWall w0;
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(-1, 0, 0), Vec3D(-getXMin(), 0, 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(1, 0, 0), Vec3D(getXMax(), 0, 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0, -1, 0), Vec3D(0, -getYMin(), 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0, 1, 0), Vec3D(0, getYMax(), 0));
        wallHandler.copyAndAddObject(w0);
        
    } else if (integerValue == 6)
    {
        PeriodicBoundary b0;
        b0.set(Vec3D(1, 0, 0), getXMin(), getXMax());
        boundaryHandler.copyAndAddObject(b0);
        b0.set(Vec3D(0, 1, 0), getYMin(), getYMax());
        boundaryHandler.copyAndAddObject(b0);
        b0.set(Vec3D(0, 0, 1), getZMin(), getZMax());
        boundaryHandler.copyAndAddObject(b0);
    } else
    {
        std::cerr << "Error in par.ini: line 1, value 1 is " << integerValue << std::endl;
        exit(-1);
    }
    
    //   2: integer (0|1) dont use | use the search pattern for linked cells
    file >> integerValue; //ignore

    //   3: real - gravity in z direction: positive points upwards
    file >> doubleValue;
    setGravity(Vec3D(0.0, 0.0, doubleValue));
    
    // line 2 =============================================================
    // Example: -10000 .5e-2
    //   1: time end of simulation - (negative resets start time to zero
    //                                and uses -time as end-time)   
    file >> doubleValue;
    if (doubleValue < 0)
        setTime(0);
    setTimeMax(fabs(doubleValue));
    
    //   2: time-step of simulation
    file >> doubleValue;
    setTimeStep(doubleValue);
    
    // line 3 =============================================================
    // Example: 1e-1 100
    file >> doubleValue;
    if (doubleValue >= 0)
    {
        //   1: time-step for output on time-series protocoll file  -> "ene"
        unsigned int savecount = static_cast<unsigned int>(round(doubleValue / getTimeStep()));
        setSaveCount(savecount);
        
        //   2: time-step for output on film (coordinate) file      -> "c3d"
        //      (fstat-output is coupled to c3d-output time-step)
        file >> doubleValue;
        savecount = static_cast<unsigned int>(round(doubleValue / getTimeStep()));
        dataFile.setSaveCount(savecount);
        fStatFile.setSaveCount(savecount);
    } else
    {
        //  or: ---------------------------------------------------------------
        //   1: negative number is multiplied to the previous log-output-time
        //   2: requires initial log-output time
        //   3: negative number is multiplied to the previous film-output-time
        //   4: requires initial film-output time
        std::cerr << "Error in par.ini: line 3, value 1 is " << doubleValue << std::endl;
        exit(-1);
    }
    
    // line 4 =============================================================
    // Example: 2000 1e5 1e3 1e2
    //   1: particle density (mass=4/3*constants::pi*density*rad^3)
    file >> doubleValue;

    //clear species handler
    speciesHandler.clear();
    auto S = new LinearViscoelasticSlidingFrictionSpecies;
    speciesHandler.addObject(S);

    setParticleDimensions(3);
    S->setDensity(doubleValue);

    //   2: linear spring constant
    file >> doubleValue;
    S->setStiffness(doubleValue);

    //   3: linear dashpot constant
    file >> doubleValue;
    S->setDissipation(doubleValue);

    //   4: background damping dashpot constant
    file >> doubleValue;
    if (doubleValue != 0.0)
        std::cerr << "Warning in par.ini: ignored background damping " << doubleValue << std::endl;
    
    // line 5 =============================================================
    // Example: 0 0
    //   1: growth rate:  d(radius) = xgrow * dt
    file >> doubleValue;
    if (doubleValue != 0.0)
        std::cerr << "Warning in par.ini: ignored growth rate " << doubleValue << std::endl;
    
    //   2: target volume_fraction
    file >> doubleValue;
    if (doubleValue != 0.0)
        std::cerr << "Warning in par.ini: ignored target volume_fraction " << doubleValue << std::endl;
    
    file.close();
    //std::cout << "Loaded par.ini file " << filename << std::endl;
    return true;
}

bool DPMBase::findNextExistingDataFile(Mdouble tMin, bool verbose)
{
    if (dataFile.getFileType() == FileType::MULTIPLE_FILES || dataFile.getFileType() == FileType::MULTIPLE_FILES_PADDED)
    {
        while (true)// This true corresponds to the if s
        {
            dataFile.openNextFile();
            //check if file exists and contains data
            int N;
            dataFile.getFstream() >> N;
            if (dataFile.getFstream().eof() || dataFile.getFstream().peek() == -1)
            {
                std::cout << "file " << dataFile.getName() << " not found" << std::endl;
                return false;
            }
            //check if tmin condition is satisfied
            Mdouble t;
            dataFile.getFstream() >> t;
            if (t > tMin)
            {
                //set_file_counter(get_file_counter()-1);
                return true;
            }
            if (verbose)
                std::cout << "Jumping file counter: " << dataFile.getCounter() << std::endl;
        }
    }
    return true;
}

bool DPMBase::readNextDataFile(unsigned int format)
{
    dataFile.openNextFile(std::fstream::in);
    //fStatFile.openNextFile();
    //Set the correct formation based of dimension if the formation is not specified by the user
    if (format == 0)
    {
        switch (getSystemDimensions())
        {
            case 1:
            case 2:
                format = 8;
                break;
            case 3:
                format = 14;
                break;
        }
        //end case

    }
    //end if

    unsigned int N = 0;
    Mdouble dummy;
    dataFile.getFstream() >> N;
    //read first parameter and check if we reached the end of the file
    if (N == 0)
    {
        //std::cout << "Reached the end of the data file" << std::endl;
        return false;
    }

    //create particles if N is not fixed
    if (particleHandler.getNumberOfObjects() != N)
    {
        if (particleHandler.getNumberOfObjects() == 0)
        {
            BaseParticle p;
            p.setSpecies(speciesHandler.getObject(0));
            particleHandler.copyAndAddObject(p);
            std::cout << "Warning in readNextDataFile: particleHandler is empty, "
                    "so newly created particles will be of BaseParticle type." << std::endl;
        }
        if (particleHandler.getNumberOfObjects() > N)
        {
            while (particleHandler.getNumberOfObjects() != N)
            {
                particleHandler.removeLastObject();
                //std::cout << "S" << particleHandler.getNumberOfObjects() << std::endl;
            }
        } else
        {
            while (particleHandler.getNumberOfObjects() != N)
            {
                particleHandler.copyAndAddObject(particleHandler.getLastObject());
                particleHandler.getLastObject()->unfix();
                //std::cout << particleHandler.getNumberOfObjects() << N << std::endl;
            }
        }
    }
    //for (auto p: particleHandler)
    //    std::cout << *p << std::endl;
#ifdef DEBUG_OUTPUT
    std::cout << "Number of particles read from file "<<N << std::endl;
#endif

    //read all other data available for the time step
    switch (format)
    {
        //This is the format_ that Stefan's and Vitaly's code saves in - note the axis rotation_
        case 7:
        {
            std::stringstream line;
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> dummy;
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> dummy;
            setTime(dummy);
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> dummy;
            setXMin(dummy);
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> dummy;
            setYMin(dummy);
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> dummy;
            setZMin(dummy);
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> dummy;
            setXMax(dummy);
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> dummy;
            setYMax(dummy);
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> dummy;
            setZMax(dummy);
            //std::cout << " time " << t <<  std::endl;
            Mdouble radius;
            Vec3D position, velocity;
            for (BaseParticle* p : particleHandler)
            {
                helpers::getLineFromStringStream(dataFile.getFstream(), line);
                line >> position.X >> position.Z >> position.Y >> velocity.X >> velocity.Z >> velocity.Y >> radius
                     >> dummy;
                p->setPosition(position);
                p->setVelocity(velocity);
                p->setOrientation(Vec3D(0.0, 0.0, 0.0));
                p->setAngularVelocity(Vec3D(0.0, 0.0, 0.0));
                p->setRadius(radius);
            }
            break;
        }
            //This is a 2D format_
        case 8:
        {
            std::stringstream line;
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> time_ >> xMin_ >> yMin_ >> xMax_ >> yMax_;
            zMin_ = 0.0;
            zMax_ = 0.0;
            Mdouble radius;
            Vec3D position, velocity, angle, angularVelocity;
            for (BaseParticle* p : particleHandler)
            {
                helpers::getLineFromStringStream(dataFile.getFstream(), line);
                line >> position.X >> position.Y >> velocity.X >> velocity.Y >> radius >> angle.Z >> angularVelocity.Z
                     >> dummy;
                p->setPosition(position);
                p->setVelocity(velocity);
                p->setOrientation(-angle);
                p->setAngularVelocity(-angularVelocity);
                p->setRadius(radius);
            } //end for all particles
            break;
        }
            //This is a 3D format_
        case 14:
        {
            std::stringstream line;
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> time_ >> xMin_ >> yMin_ >> zMin_ >> xMax_ >> yMax_ >> zMax_;
            Mdouble radius;
            Vec3D position, velocity, angle, angularVelocity;
            for (BaseParticle* p : particleHandler)
            {
                helpers::getLineFromStringStream(dataFile.getFstream(), line);
                line >> position >> velocity >> radius >> angle >> angularVelocity;
                p->setPosition(position);
                p->setVelocity(velocity);
                p->setOrientation(angle);
                p->setAngularVelocity(angularVelocity);
                p->setRadius(radius);
            } //end read into existing particles
            logger(INFO, "read % particles", particleHandler.getNumberOfObjects());
            break;
        }
            //This is a 3D format_
        case 15:
        {
            std::stringstream line;
            helpers::getLineFromStringStream(dataFile.getFstream(), line);
            line >> time_ >> xMin_ >> yMin_ >> zMin_ >> xMax_ >> yMax_ >> zMax_;
            Mdouble radius;
            Vec3D position, velocity, angle, angularVelocity;
            for (BaseParticle* p : particleHandler)
            {
                helpers::getLineFromStringStream(dataFile.getFstream(), line);
                line >> position >> velocity >> radius >> angle >> angularVelocity >> dummy >> dummy;
                p->setPosition(position);
                p->setVelocity(velocity);
                p->setOrientation(angle);
                p->setAngularVelocity(angularVelocity);
                p->setRadius(radius);
            } //end for all particles
            break;
        }
    } //end switch statement
    particleHandler.computeAllMasses();
    return true;
}

void DPMBase::writeRestartFile()
{
    write(restartFile.getFstream());
}

bool DPMBase::readRestartFile(bool restarted)
{
    if (restartFile.open(std::fstream::in))
    {
        read(restartFile.getFstream());
        restartFile.close();
        setRestarted(restarted);
        return true;
    } else {
        //std::cerr << restartFile.getName() << " could not be loaded" << std::endl;
        return false;
    }
}

int DPMBase::readRestartFile(std::string fileName)
{
    restartFile.setName(fileName);
    return readRestartFile();
}

void DPMBase::computeInternalForces(BaseParticle* P1, BaseParticle* P2)
{
    //this is because the rough bottom allows overlapping fixed particles
    if (P1->isFixed() && P2->isFixed())
    {
        return;
    }

    if (P1->getPeriodicFromParticle() && P2->getPeriodicFromParticle())
    {
        return;
    }

    BaseParticle* PI, * PJ;
    if (P1->getId() > P2->getId())
    {
        PI = P2;
        PJ = P1;
    } else
    {
        PI = P1;
        PJ = P2;
    }

    //if statement above ensures that the PI has the lower id than PJ
    std::vector<BaseInteraction*> interactions = PJ->getInteractionWith(PI, getTime() + getTimeStep(),
                                                                        &interactionHandler);
    for (auto i : interactions)
    {
        i->computeForce();

        PI->addForce(i->getForce());
        PJ->addForce(-i->getForce());

        if (getRotation())
        {
            PI->addTorque(i->getTorque() - Vec3D::cross(PI->getPosition() - i->getContactPoint(), i->getForce()));
            PJ->addTorque(-i->getTorque() + Vec3D::cross(PJ->getPosition() - i->getContactPoint(), i->getForce()));
        }
    }
}

void DPMBase::computeExternalForces(BaseParticle* CI)
{
    if (!CI->isFixed())
    {
        // Gravitational force
        CI->addForce(getGravity() * CI->getMass());
        // Still calls this in compute External Forces.
        computeForcesDueToWalls(CI);
    }
}

void DPMBase::computeForcesDueToWalls(BaseParticle* pI)
{
    //No need to compute interactions between periodic particle images and walls
    if (pI->getPeriodicFromParticle() != nullptr)
        return;

    for (BaseWall* w : wallHandler)
    {
        std::vector<BaseInteraction*> interactions = w->getInteractionWith(pI, getTime() + getTimeStep(),
                                                                           &interactionHandler);

        for (auto i : interactions)
        {
            i->computeForce();

            pI->addForce(i->getForce());
            w->addForce(-i->getForce());

            if (getRotation()) // getRotation() returns a boolean.
            {
                pI->addTorque(i->getTorque() - Vec3D::cross(pI->getPosition() - i->getContactPoint(), i->getForce()));
                w->addTorque(-i->getTorque() + Vec3D::cross(w->getPosition() - i->getContactPoint(), i->getForce()));
            }
        }
    }
}

void DPMBase::integrateBeforeForceComputation()
{
    for_each(particleHandler.begin(), particleHandler.end(), [this](BaseParticle* p)
    {
        p->integrateBeforeForceComputation(getTime(), getTimeStep());
    });
    for_each(wallHandler.begin(), wallHandler.end(), [this](BaseWall* w)
    {
        w->integrateBeforeForceComputation(getTime(), getTimeStep());
    });
}

void DPMBase::checkInteractionWithBoundaries()
{
    for (BaseBoundary* b : boundaryHandler)
    {
        for (auto p = particleHandler.begin(); p != particleHandler.end(); ++p)
        {
            if (b->checkBoundaryAfterParticleMoved(*p, particleHandler))
                //Returns true if the particle is deleted
                --p;
        }
    }
}

void DPMBase::integrateAfterForceComputation()
{
    for_each(particleHandler.begin(), particleHandler.end(), [this](BaseParticle* p)
    {
        p->integrateAfterForceComputation(getTime(), getTimeStep());
    });
    for_each(wallHandler.begin(), wallHandler.end(), [this](BaseWall* w)
    {
        w->integrateAfterForceComputation(getTime(), getTimeStep());
    });
}

//void DPMBase::statisticsFromRestartData(const char *name)
//{
//    ///todo{Check this whole function}
//    //This function loads all MD data
//    readRestartFile();
//
//    //This creates the file statistics will be saved to
//    std::stringstream ss("");
//    ss << name << ".stat";
//    statFile.setName(ss.str());
//    statFile.setOpenMode(std::fstream::out);
//    statFile.open();
//
//    // Sets up the initial conditions for the simulation
//    // setupInitialConditions();
//    // Setup the previous position arrays and mass of each particle.
//    computeParticleMasses();
//    // Other routines required to jump-start the simulation
//    actionsBeforeTimeLoop();
//    initialiseStatistics();
//    hGridActionsBeforeTimeLoop();
//    writeEneHeader(eneFile.getFstream());
//
//    while (readDataFile(dataFile.getName().c_str()))
//    {
//        hGridActionsBeforeTimeLoop();
//        actionsBeforeTimeStep();
//        checkAndDuplicatePeriodicParticles();
//        hGridActionsBeforeTimeStep();
//        computeAllForces();
//        removeDuplicatePeriodicParticles();
//        actionsAfterTimeStep();
//        writeEneTimestep(eneFile.getFstream());
//        std::cout << std::setprecision(6) << getTime() << std::endl;
//    }
//
//    dataFile.close();
//    statFile.close();
//}

void DPMBase::computeAllForces()
{
    for (BaseParticle* const p : particleHandler)
    {
        p->setForce(Vec3D(0.0, 0.0, 0.0));
        p->setTorque(Vec3D(0.0, 0.0, 0.0));
    }
    for (BaseWall* const w : wallHandler)
    {
        w->setForce(Vec3D(0.0, 0.0, 0.0));
        w->setTorque(Vec3D(0.0, 0.0, 0.0));
    } //end reset forces loop

#ifdef DEBUG_OUTPUT
    std::cerr << "Have all forces set to zero " << std::endl;
#endif

    
    for (BaseParticle* p : particleHandler)
    {

        computeInternalForces(p);
        //end inner loop over contacts.

        computeExternalForces(p);

    }

#ifdef CONTACT_LIST_HGRID
    //possibleContactList.write(std::cout);
    PossibleContact* Curr=possibleContactList.getFirstPossibleContact();
    while(Curr)
    {   
        computeInternalForces(Curr->getP1(),Curr->getP2());
        Curr=Curr->getNext();
    }
#endif
    //end outer loop over contacts.

}

void DPMBase::computeInternalForces(BaseParticle* i)
{
    broadPhase(i);
}

void DPMBase::write(std::ostream& os, bool writeAllParticles) const
{
    os << "restart_version " << "1.0";
    FilesAndRunNumber::write(os);
    os << "xMin " << getXMin()
       << " xMax " << getXMax()
       << " yMin " << getYMin()
       << " yMax " << getYMax()
       << " zMin " << getZMin()
       << " zMax " << getZMax() << std::endl
       << "timeStep " << getTimeStep()
       << " time " << getTime()
       << " ntimeSteps " << ntimeSteps_
       << " timeMax " << getTimeMax() << std::endl
       << "systemDimensions " << getSystemDimensions()
       << " particleDimensions " << getParticleDimensions()
       << " gravity " << getGravity();
    os << " writeVTK " << writeParticlesVTK_ << " " << writeWallsVTK_ << " " << interactionHandler.getWriteVTK();
    os << " random "; random.write(os);
    //only write xBallsArguments if they are nonzero
    if (getXBallsAdditionalArguments().compare(""))
        os << " xBallsArguments " << getXBallsAdditionalArguments();
    os << std::endl;

    speciesHandler.write(os);
    os << "Walls " << wallHandler.getNumberOfObjects() << std::endl;
    for (BaseWall* w : wallHandler)
        os << (*w) << std::endl;
    os << "Boundaries " << boundaryHandler.getNumberOfObjects() << std::endl;
    for (BaseBoundary* b : boundaryHandler)
        os << (*b) << std::endl;
    if (writeAllParticles || particleHandler.getNumberOfObjects() < 4)
    {
        particleHandler.write(os);
    } else
    {
        os << "Particles " << particleHandler.getNumberOfObjects() << std::endl;
        for (unsigned int i = 0; i < 2; i++)
            os << *particleHandler.getObject(i) << std::endl;
        os << "..." << std::endl;
    }
    if (writeAllParticles || interactionHandler.getNumberOfObjects() < 4)
    {
        interactionHandler.write(os);
    } else
    {
        os << "Interactions " << interactionHandler.getNumberOfObjects() << std::endl;
        for (unsigned int i = 0; i < 2; i++)
            os << *interactionHandler.getObject(i) << std::endl;
        os << "..." << std::endl;
    }
}

void DPMBase::read(std::istream& is)
{
    std::string dummy;
    is >> dummy;
    if (strcmp(dummy.c_str(), "restart_version"))
    {
        //only very old files did not have a restart_version
        logger.log(Log::FATAL, "Error in DPMBase::read(is): this is not a valid restart file");
    } else
    {
        is >> restartVersion_;
        if (!restartVersion_.compare("1.0"))
        {
            FilesAndRunNumber::read(is);
            is >> dummy >> xMin_
               >> dummy >> xMax_
               >> dummy >> yMin_
               >> dummy >> yMax_
               >> dummy >> zMin_
               >> dummy >> zMax_
               >> dummy >> timeStep_
               >> dummy >> time_
               >> dummy >> ntimeSteps_
               >> dummy >> timeMax_
               >> dummy >> systemDimensions_
               >> dummy >> particleDimensions_
               >> dummy >> gravity_;

            std::stringstream line(std::stringstream::in | std::stringstream::out);
            helpers::getLineFromStringStream(is, line);
            line >> dummy;
            if (!dummy.compare("writeVTK")) {
                FileType writeInteractionsVTK = FileType::NO_FILE;
                line >> writeParticlesVTK_ >> writeWallsVTK_ >> writeInteractionsVTK >> dummy;
                interactionHandler.setWriteVTK(writeInteractionsVTK);
            }
            if (!dummy.compare("random"))
            {
                helpers::getLineFromStringStream(line, line);
                random.read(line);
                line >> dummy;
                helpers::getLineFromStringStream(line, line);
                setXBallsAdditionalArguments(line.str());
            }

            speciesHandler.read(is);

            unsigned int N;
            //std::stringstream line(std::stringstream::in | std::stringstream::out);
            is >> dummy >> N;
            if (dummy.compare("Walls"))
                logger(ERROR, "DPMBase::read(is): Error during restart: 'Walls' argument could not be found.");
            wallHandler.clear();
            wallHandler.setStorageCapacity(N);
            helpers::getLineFromStringStream(is, line);

            for (unsigned int i = 0; i < N; i++)
            {
                helpers::getLineFromStringStream(is, line);
                wallHandler.readObject(line);
            }

            is >> dummy >> N;
            boundaryHandler.clear();
            boundaryHandler.setStorageCapacity(N);
            if (dummy.compare("Boundaries"))
                logger(ERROR, "DPMBase::read(is): Error during restart: 'Boundaries' argument could not be found.");
            helpers::getLineFromStringStream(is, line);
            for (unsigned int i = 0; i < N; i++)
            {
                helpers::getLineFromStringStream(is, line);
                boundaryHandler.readObject(line);
            }

            is >> dummy >> N;
            if (dummy.compare("Particles"))
                logger(ERROR, "DPMBase::read(is): Error during restart: 'Boundaries' argument could not be found.");
            particleHandler.clear();
            particleHandler.setStorageCapacity(N);
            helpers::getLineFromStringStream(is, line);
            for (unsigned int i = 0; i < N; i++)
            {
                helpers::getLineFromStringStream(is, line);
                particleHandler.readObject(line);
                //particleHandler.getLastObject()->computeMass();
            }

            interactionHandler.read(is);

        } else if (!restartVersion_.compare("3"))
        {
            logger.log(Log::INFO, "DPMBase::read(is): restarting from an old restart file (restart_version %).",
                       restartVersion_);
            readOld(is);
        } else
        {
            //only very old files did not have a restart_version
            logger.log(Log::FATAL,
                       "Error in DPMBase::read(is): restart_version % cannot be read; use an older version of Mercury to upgrade the file",
                       restartVersion_);
        }
    }
}

void DPMBase::readOld(std::istream& is)
{
    std::string dummy;
    is >> dummy >> dummy;
    setName(dummy);

    unsigned int saveCountData, saveCountEne, saveCountStat, saveCountFStat;
    unsigned int fileTypeFstat, fileTypeData, fileTypeEne, fileTypeRestart;
    is >> dummy >> xMin_
       >> dummy >> xMax_
       >> dummy >> yMin_
       >> dummy >> yMax_
       >> dummy >> zMin_
       >> dummy >> zMax_
       >> dummy >> timeStep_
       >> dummy >> time_
       >> dummy >> timeMax_
       >> dummy >> saveCountData
       >> dummy >> saveCountEne
       >> dummy >> saveCountStat
       >> dummy >> saveCountFStat
       >> dummy >> systemDimensions_
       >> dummy >> gravity_
       >> dummy >> fileTypeFstat
       >> dummy >> fileTypeData
       >> dummy >> fileTypeEne;
    dataFile.setSaveCount(saveCountData);
    eneFile.setSaveCount(saveCountEne);
    statFile.setSaveCount(saveCountStat);
    fStatFile.setSaveCount(saveCountFStat);
    
    fStatFile.setFileType(static_cast<FileType>(fileTypeFstat));
    dataFile.setFileType(static_cast<FileType>(fileTypeData));
    eneFile.setFileType(static_cast<FileType>(fileTypeEne));
    
    //this is optional to allow restart files with and without restartFile.getFileType()
    is >> dummy;
    if (!strcmp(dummy.c_str(), "options_restart"))
    {
        is >> fileTypeRestart;
        restartFile.setFileType(static_cast<FileType>(fileTypeRestart));
    }

    speciesHandler.read(is);
    wallHandler.read(is);
    boundaryHandler.read(is);
    particleHandler.read(is);
}

void DPMBase::checkSettings()
{
    if (speciesHandler.getNumberOfObjects() == 0)
    {
        std::cerr << "Error in solve(): at least one species has to be defined" << std::endl;
        exit(-1);
    }
    if (getParticleDimensions() == 0)
    {
        std::cerr << "Error in solve(): particleDimensions not specified" << std::endl;
        exit(-1);
    }
    if (getSystemDimensions() == 0)
    {
        std::cerr << "Error in solve(): systemDimensions not specified" << std::endl;
        exit(-1);
    }
    if (getTimeStep() == 0.0)
    {
        std::cerr << "Error in solve(): timeStep not specified" << std::endl;
        exit(-1);
    }
}

void DPMBase::writeOutputFiles()
{
    if (fStatFile.saveCurrentTimestep(ntimeSteps_))
        writeFstatHeader(fStatFile.getFstream());

    if (eneFile.saveCurrentTimestep(ntimeSteps_))
    {
        if (eneFile.getCounter() == 1 || eneFile.getFileType() == FileType::MULTIPLE_FILES ||
            eneFile.getFileType() == FileType::MULTIPLE_FILES_PADDED)
            writeEneHeader(eneFile.getFstream());
        writeEneTimestep(eneFile.getFstream());
    }

    if (dataFile.saveCurrentTimestep(ntimeSteps_))
    {
        printTime();
        if ((getRestarted() || dataFile.getCounter() == 1) && dataFile.getFileType() != FileType::NO_FILE)
            writeXBallsScript();
        outputXBallsData(dataFile.getFstream());
        writeVTK();
    }

//    if (statFile.saveCurrentTimestep(ntimeSteps_))
//    {
//        outputStatistics();
//        gatherContactStatistics();
//        processStatistics(true);
//    }

    //write restart file last, otherwise the output cunters are wrong
    if (restartFile.saveCurrentTimestep(ntimeSteps_))
    {
        writeRestartFile();
        restartFile.close(); //overwrite old restart file if FileType::ONE_FILE
    }
}

void DPMBase::solve()
{
    logger(DEBUG, "Entered solve");
#ifdef CONTACT_LIST_HGRID
    logger(INFO,"Using CONTACT_LIST_HGRID");
#endif
    
    if (!getRestarted())
    {
        ntimeSteps_ = 0;
        resetFileCounter();
        setTime(0.0);
        setNextSavedTimeStep(0); //reset the counter
        //this is to ensure that the interaction time stamps agree with the resetting of the time value
        for (auto& i : interactionHandler)
            i->setTimeStamp(0);
        setupInitialConditions();
        logger(DEBUG, "Have created the particles initial conditions");
    }
    else
    {
        actionsOnRestart();
    }

    checkSettings();

    if (getRunNumber() > 0 && !getRestarted())
    {
        std::stringstream name;
        name << getName() << "." << getRunNumber();
        setName(name.str());
    }

    //append_ determines if files have to be appended or not
    if (getAppend())
    {
        setOpenMode(std::fstream::out | std::fstream::app);
        restartFile.setOpenMode(std::fstream::out); //Restart files should always be overwritten.
    }
    else
    {
        setOpenMode(std::fstream::out);
    }

    initialiseStatistics();

    // Setup the mass of each particle.
    particleHandler.computeAllMasses();

    // Other initialisations
    //max_radius = getLargestParticle()->getRadius();
    actionsBeforeTimeLoop();
    hGridActionsBeforeTimeLoop();

    // do a first force computation
    checkAndDuplicatePeriodicParticles();
    hGridActionsBeforeTimeStep();
    time_ -= timeStep_; //to enforce that old  time stamps are created
    ntimeSteps_--;
    computeAllForces();
    ntimeSteps_++;
    time_ += timeStep_;
    removeDuplicatePeriodicParticles();
    interactionHandler.actionsAfterTimeStep();
    logger(DEBUG, "Have computed the initial values for the forces ");

    // This is the main loop over advancing time
    while (getTime() < getTimeMax() && continueSolve())
    {
        writeOutputFiles(); //everything is written at the beginning of the timestep!
        // Loop over all particles doing the time integration step
        hGridActionsBeforeIntegration();
        integrateBeforeForceComputation();
        checkInteractionWithBoundaries(); // INSERTION boundaries handled
        hGridActionsAfterIntegration();

        // Compute forces

        // INSERTION/DELETION boundary flag change
        for (BaseBoundary* b : boundaryHandler)
        {
            b->checkBoundaryBeforeTimeStep(this);
        }

        actionsBeforeTimeStep();

        checkAndDuplicatePeriodicParticles();

        hGridActionsBeforeTimeStep();

        computeAllForces();

        removeDuplicatePeriodicParticles();

        actionsAfterTimeStep();

        // Loop over all particles doing the time integration step
        hGridActionsBeforeIntegration();
        integrateAfterForceComputation();

        checkInteractionWithBoundaries(); // DELETION boundaries handled
        hGridActionsAfterIntegration();

        //erase interactions that have not been used during the last timestep
        interactionHandler.eraseOldInteractions(getTime());
        interactionHandler.actionsAfterTimeStep();

        time_ += timeStep_;
        ntimeSteps_++;
    }
    //force writing of the last time step
    setNextSavedTimeStep(ntimeSteps_);
    writeOutputFiles();

    //end loop over interaction count
    actionsAfterSolve();

    //To make sure getTime gets the correct time for outputting statistics
    finishStatistics();

    closeFiles();
}

bool DPMBase::readArguments(int argc, char* argv[])
{
    bool isRead = true;
    for (int i = 1; i < argc; i += 2)
    {
        std::cout << "interpreting input argument " << argv[i];
        for (int j = i + 1; j < argc; j++)
        {
            if (argv[j][0] == '-')
                break;
            std::cout << " " << argv[j];
        }
        std::cout << std::endl;
        isRead &= readNextArgument(i, argc, argv);
        if (!isRead)
        {
            std::cerr << "Warning: not all arguments read correctly!" << std::endl;
            exit(-1);
        }
    }
    return isRead;
}

bool DPMBase::readNextArgument(int& i, int argc, char* argv[])
{
    if (!strcmp(argv[i], "-name"))
    {
        setName(argv[i + 1]);
    } else if (!strcmp(argv[i], "-xmin"))
    {
        setXMin(atof(argv[i + 1]));
    } else if (!strcmp(argv[i], "-ymin"))
    {
        setYMin(atof(argv[i + 1]));
    } else if (!strcmp(argv[i], "-zmin"))
    {
        setZMin(atof(argv[i + 1]));
    } else if (!strcmp(argv[i], "-xmax"))
    {
        setXMax(atof(argv[i + 1]));
    } else if (!strcmp(argv[i], "-ymax"))
    {
        setYMax(atof(argv[i + 1]));
    } else if (!strcmp(argv[i], "-zmax"))
    {
        setZMax(atof(argv[i + 1]));
        //} else if (!strcmp(argv[i],"-svn")) {
        //  std::cout << "svn version " << SVN_VERSION << std::endl;
        //  i--;
    } else if (!strcmp(argv[i], "-addS"))
    {
        setAddSmall(atof(argv[i + 1]));
    } else if (!strcmp(argv[i], "-addL"))
    {
        setAddLarge(atof(argv[i + 1]));
    } else if (!strcmp(argv[i], "-dt"))
    {
        Mdouble old = getTimeStep();
        setTimeStep(atof(argv[i + 1]));
        std::cout << "  reset dt from " << old << " to " << getTimeStep() << std::endl;
    }
//    else if (!strcmp(argv[i], "-Hertz"))
//    {
//        speciesHandler.getObject(0)->setForceType(ForceType::HERTZ);
//        i--;
//    }
    else if (!strcmp(argv[i], "-tmax"))
    {
        Mdouble old = getTimeMax();
        setTimeMax(atof(argv[i + 1]));
        std::cout << "  reset timeMax from " << old << " to " << getTimeMax() << std::endl;
    } else if (!strcmp(argv[i], "-saveCount"))
    {
        Mdouble old = dataFile.getSaveCount();
        setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
        setNextSavedTimeStep(getNtimeSteps() + restartFile.getSaveCount());
        std::cout << "  reset saveCount from " << old << " to " << dataFile.getSaveCount() << std::endl;
    } else if (!strcmp(argv[i], "-saveCountData"))
    {
        dataFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-saveCountFStat"))
    {
        fStatFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-saveCountStat"))
    {
        statFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-saveCountEne"))
    {
        eneFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-saveCountRestart"))
    {
        restartFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-dim"))
    {
        setSystemDimensions(static_cast<unsigned int>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-gravity"))
    {
        setGravity(Vec3D(atof(argv[i + 1]), atof(argv[i + 2]), atof(argv[i + 3])));
        i += 2;
    } else if (!strcmp(argv[i], "-fileType"))
    { //uses int input
        setFileType(static_cast<FileType>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-fileTypeFStat"))
    { //uses int input
        fStatFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-fileTypeRestart"))
    {
        restartFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-fileTypeData"))
    {
        dataFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-fileTypeStat"))
    {
        statFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-fileTypeEne"))
    {
        eneFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
    } else if (!strcmp(argv[i], "-auto_number"))
    {
        autoNumber();
        i--;
    }
//    else if (!strcmp(argv[i], "-number_of_saves"))
//    {
//        set_number_of_saves_all(atof(argv[i + 1]));
//    }
    else if (!strcmp(argv[i], "-restart") || !strcmp(argv[i], "-r"))
    {
        std::string filename;
        
        //use default filename if no argument is given
        if (i + 1 >= argc || argv[i + 1][0] == '-')
        {
            i--;
            filename = getName();
            std::cout << getName() << std::endl;
        } else
            filename = argv[i + 1];

        //add ".restart" if necessary
        if (filename.find(".restart") == std::string::npos)
        {
            filename = filename + ".restart";
        }
        
        std::cout << "restart from " << filename << std::endl;
        readRestartFile(filename);
    } else if (!strcmp(argv[i], "-clean") || !strcmp(argv[i], "-c"))
    {
        std::cout << "Remove old " << getName() << ".* files" << std::endl;
//        std::string filename;
        std::ostringstream filename;
        std::vector<std::string> ext{".restart", ".stat", ".fstat", ".data", ".ene", ".xballs"};
        for (unsigned int j = 0; j < ext.size(); j++)
        {
            // remove files with given extension for FileType::ONE_FILE
            filename.clear();
            filename << getName() << ext[j];
            if (!remove(filename.str().c_str()))
            {
                std::cout << "  File " << filename.str() << " successfully deleted" << std::endl;
            }
            // remove files with given extension for FileType::MULTIPLE_FILES
            unsigned int k = 0;
            filename.clear();
            filename << getName() << ext[j] << '.' << k;
            while (!remove(filename.str().c_str()))
            {
                std::cout << "  File " << filename.str() << " successfully deleted" << std::endl;
                filename.clear();
                filename << getName() << ext[j] << '.' << ++k;
            }
            // remove files with given extension for FileType::MULTIPLE_FILES_PADDED
            k = 0;
            filename.clear();
            filename << getName() << ext[j] << '.' << to_string_padded(k);
            while (!remove(filename.str().c_str()))
            {
                std::cout << "  File " << filename.str() << " successfully deleted" << std::endl;
                filename.clear();
                filename << getName() << ext[j] << '.' << to_string_padded(++k);
            }
        }
        i--;
    } else if (!strcmp(argv[i], "-data"))
    {
        std::string filename = argv[i + 1];
        readDataFile(filename.c_str());
    }
//    else if (!strcmp(argv[i], "-k"))
//    {
//        Mdouble old = getSpecies(0)->getStiffness();
//        getSpecies(0)->setStiffness(atof(argv[i + 1]));
//        std::cout << "  reset k from " << old << " to " << getSpecies(0)->getStiffness() << std::endl;
//    }
//    else if (!strcmp(argv[i], "-dissipation") || !strcmp(argv[i], "-disp"))
//    {
//        Mdouble old = getSpecies(0)->getDissipation();
//        getSpecies(0)->setDissipation(atof(argv[i + 1]));
//        std::cout << "  reset getDissipation() from " << old << " to " << getSpecies(0)->getDissipation() << std::endl;
//    }
//    else if (!strcmp(argv[i], "-kt"))
//    {
//        Mdouble old = getSpecies(0)->getSlidingStiffness();
//        getSpecies(0)->setSlidingStiffness(atof(argv[i + 1]));
//        std::cout << "  reset kt from " << old << " to " << getSpecies(0)->getSlidingStiffness() << std::endl;
//    }
//    else if (!strcmp(argv[i], "-dispt"))
//    {
//        Mdouble old = getSpecies(0)->getSlidingDissipation();
//        getSpecies(0)->setSlidingDissipation(atof(argv[i + 1]));
//        std::cout << "  reset dispt from " << old << " to " << getSpecies(0)->getSlidingDissipation() << std::endl;
//    }
//    else if (!strcmp(argv[i], "-krolling"))
//    {
//        Mdouble old = getSpecies(0)->getRollingStiffness();
//        getSpecies(0)->setRollingStiffness(atof(argv[i + 1]));
//        std::cout << "  reset krolling from " << old << " to " << getSpecies(0)->getRollingStiffness() << std::endl;
//    }
//    else if (!strcmp(argv[i], "-disprolling"))
//    {
//        Mdouble old = getSpecies(0)->getRollingDissipation();
//        getSpecies(0)->setRollingDissipation(atof(argv[i + 1]));
//        std::cout << "  reset disprolling from " << old << " to " << getSpecies(0)->getRollingDissipation() << std::endl;
//    }
//    else if (!strcmp(argv[i], "-mu"))
//    {
//        Mdouble old = getSpecies(0)->getSlidingFrictionCoefficient();
//        getSpecies(0)->setSlidingFrictionCoefficient(atof(argv[i + 1]));
//        std::cout << "  reset mu from " << old << " to " << getSpecies(0)->getSlidingFrictionCoefficient() << std::endl;
//    }
//    else if (!strcmp(argv[i], "-murolling"))
//    {
//        Mdouble old = getSpecies(0)->getRollingFrictionCoefficient();
//        getSpecies(0)->setRollingFrictionCoefficient(atof(argv[i + 1]));
//        std::cout << "  reset murolling from " << old << " to " << getSpecies(0)->getRollingFrictionCoefficient() << std::endl;
//    }
    else if (!strcmp(argv[i], "-randomise") || !strcmp(argv[i], "-randomize"))
    {
        random.randomise();
        i--;
    }
//    else if (!strcmp(argv[i], "-k0"))
//    {
//        Mdouble old = speciesHandler.getObject(0)->getAdhesionStiffness();
//        speciesHandler.getObject(0)->setAdhesionStiffness(atof(argv[i + 1]));
//        std::cout << "  reset k0 from " << old << " to " << speciesHandler.getObject(0)->getAdhesionStiffness() << std::endl;
//    }
//    else if (!strcmp(argv[i], "-f0"))
//    {
//        Mdouble old = speciesHandler.getObject(0)->getBondForceMax();
//        speciesHandler.getObject(0)->setBondForceMax(atof(argv[i + 1]));
//        std::cout << "  reset f0 from " << old << " to " << speciesHandler.getObject(0)->getBondForceMax() << std::endl;
//    }
//    else if (!strcmp(argv[i], "-AdhesionForceType"))
//    {
//        AdhesionForceType old = speciesHandler.getObject(0)->getAdhesionForceType();
//        speciesHandler.getObject(0)->setAdhesionForceType(argv[i + 1]);
//        std::cout << "  reset AdhesionForceType from "
//                << static_cast<signed char>(old) << " to "
//                << static_cast<signed char>(speciesHandler.getObject(0)->getAdhesionForceType()) << std::endl;
//    }
    else if (!strcmp(argv[i], "-append"))
    {
        setAppend(true);
        i--;
    } else if (!strcmp(argv[i], "-fixedParticles"))
    {
        setFixedParticles(static_cast<unsigned int>(atoi(argv[i + 1])));
    }
//    else if (!strcmp(argv[i], "-rho"))
//    {
//        Mdouble old = speciesHandler.getObject(0)->getDensity();
//        speciesHandler.getObject(0)->setDensity(atof(argv[i + 1]));
//        std::cout << "  reset rho from " << old << " to " << speciesHandler.getObject(0)->getDensity() << std::endl;
//    }
//    else if (!strcmp(argv[i], "-dim_particle"))
//    {
//        setParticleDimensions(atoi(argv[i + 1]));
//    }
    else if (!strcmp(argv[i], "-counter"))
    {
        setRunNumber(atoi(argv[i + 1]));
    } else
        return false;
    return true; //returns true if argv is found
}

bool DPMBase::checkParticleForInteraction(const BaseParticle& p)
{
#ifdef MERCURY_USE_MPI
    int localInteraction = checkParticleForInteractionLocal(p);

    MPIContainer& communicator = MPIContainer::Instance();
    int numberOfProcessors = communicator.getNumberOfProcessors();

    //The root gathers all values and computes the global value
    int *interactionList = nullptr;
    if (communicator.getProcessorID() == 0)
    {
        interactionList = new int [numberOfProcessors];
    }

    //Gather all local values
    communicator.gather(localInteraction,interactionList);

    //Compute the global value
    int globalInteraction = 1;
    if (communicator.getProcessorID() == 0)
    {
        for (int i = 0; i<numberOfProcessors; i++)
        {
            if (interactionList[i] == 0)
            {
                globalInteraction = 0;
                break;
            }
        }
    }
    //The root now tells the other processors what the global value for the interaction is
    communicator.broadcast(globalInteraction);

    //Convert the result back to bool
    bool interaction = globalInteraction;
#else
    bool interaction = checkParticleForInteractionLocalPeriodic(p);
#endif
    return interaction;
}

bool DPMBase::checkParticleForInteractionLocalPeriodic(const BaseParticle& p) {
    //A vector of ghost particles of the particle that is to be inserted (empty if no periodic boundaries are present)
    std::vector<BaseParticle> pPeriodic;
    for (BaseBoundary* b : boundaryHandler) {
        PeriodicBoundary* pb = dynamic_cast<PeriodicBoundary*>(b);
        if (pb) {
            for (int i=pPeriodic.size()-1; i>=0; --i) {
                if (particleHandler.getNumberOfObjects()>0 && pb->getDistance(pPeriodic[i]) < pPeriodic[i].getInteractionRadius() + particleHandler.getLargestParticle()->getInteractionRadius()) {
                    pPeriodic.push_back(pPeriodic[i]);
                    pb->shiftPosition(&pPeriodic.back());
                }
            }
            if (particleHandler.getNumberOfObjects()>0 && pb->getDistance(p) < p.getInteractionRadius() + particleHandler.getLargestParticle()->getInteractionRadius()) {
                pPeriodic.push_back(p);
                pb->shiftPosition(&pPeriodic.back());
            }
        }
    }
    //check the particle AND the ghost particles for intersection problems
    bool insertable = checkParticleForInteractionLocal(p);
    for (BaseParticle& pP : pPeriodic) {
        insertable &= checkParticleForInteractionLocal(pP);
    }
    return insertable;
}

bool DPMBase::checkParticleForInteractionLocal(const BaseParticle& p)
{
    Mdouble distance;
    Vec3D normal;

    //Check if it has no collision with walls
    for (BaseWall* w : wallHandler)
    {
        //returns false if the function getDistanceAndNormal returns true,
        //i.e. if there exists an interaction between wall and particle
        //\todo TW getDistanceAndNormal(p,distance,normal) should ideally be replaced by a inContact(p) function, as it doesn't require distance and normal for anything (and walls now can have multiple contacts, soon particles can have it too.)
        if (w->getDistanceAndNormal(p, distance, normal))
        {
            //std::cout<<"failure: Collision with wall: "<<**it<<std::endl;
            return false;
        }
        else
        {
            //std::cout<<"No collision with wall: "<<**it<<std::endl;
        }
    }

    //Check if it has no collision with other particles
    for (BaseParticle* q : particleHandler)
    {
        //returns false if the particle separation is less than the relevant sum of interaction radii
        //(i.e. another particle is in contact with P)
        if (Vec3D::getDistanceSquared(q->getPosition(), p.getPosition())
            < mathsFunc::square(q->getInteractionRadius() + p.getInteractionRadius()))
        {
            //std::cout<<"failure: Collision with particle "<<**it<<std::endl;
            return false;
        }
        else
        {
            //std::cout<<"No collision with particle "<<**it<<std::endl;
        }
    }
    return true;
}

void DPMBase::removeDuplicatePeriodicParticles()
{
    for (unsigned int i = particleHandler.getNumberOfObjects();
         i >= 1 && particleHandler.getObject(i - 1)->getPeriodicFromParticle() != nullptr; i--)
    {
        while (particleHandler.getObject(i - 1)->getInteractions().size() > 0)
        {
            interactionHandler.removeObjectKeepingPeriodics(
                    particleHandler.getObject(i - 1)->getInteractions().front()->getIndex());
        }
        particleHandler.removeObject(i - 1);
    }
}

void DPMBase::checkAndDuplicatePeriodicParticles()
{
    for (BaseBoundary* b : boundaryHandler) // For all pointers of type BaseBoundary in Boundary handler
    {
        unsigned int N = particleHandler.getNumberOfObjects(); //Required because the number of particles increases
        for (unsigned int i = 0; i < N; i++)
        {
            b->createPeriodicParticles(particleHandler.getObject(i), particleHandler);
        }
    }
}

void DPMBase::outputInteractionDetails() const
{
    std::cout << "Interactions currently in the handler:" << std::endl;
    for (BaseInteraction* p : interactionHandler)
    {
        p->write(std::cout);
        std::cout << std::endl;
        std::cout << "Interaction " << p->getName() << " " << p->getId() << " between " << p->getP()->getId() << " and "
                  << p->getI()->getId() << std::endl;
    }
}

bool DPMBase::isTimeEqualTo(Mdouble time) const
{
    return getTime() <= time && getTime() + getTimeStep() > time;
}

Mdouble DPMBase::domainSize() const
{
    return (getXMax() - getXMin()) * (getYMax() - getYMin()) * (getZMax() - getZMin());
}