Helpers.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 "Helpers.h"
#include <cmath>
#include <sys/stat.h>
#include <sstream>
#include <string>
#include <Logger.h>
#include <DPMBase.h>
#include <Walls/InfiniteWall.h>
#include <Particles/BaseParticle.h>
#include <Species/BaseSpecies.h>
#include <assert.h>
#include "Math/ExtendedMath.h"


helpers::KAndDisp helpers::computeKAndDispFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble mass)
{
    helpers::KAndDisp ans;
    ans.k = computeKFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(tc, r, mass);
    ans.disp = computeDispFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(tc, r, mass);
    return ans;
}

Mdouble helpers::computeCollisionTimeFromKAndDispAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass)
{
    if (k <= 0)
    {
        std::cerr << "Error in helpers::computeCollisionTimeFromKandDispAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass) stiffness is not set or has an unexpected value, (stiffness=" << k << ")" << std::endl;
        exit(-1);
    }
    if (disp < 0)
    {
        std::cerr << "Error in helpers::computeCollisionTimeFromKandDispAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass) dissipation is not set or has an unexpected value, (dissipation=" << disp << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeCollisionTimeFromKandDispAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    if (4.0 * k / mass - mathsFunc::square(disp / mass) < 0)
    {
        std::cerr << "Error in helpers::computeCollisionTimeFromKandDispAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass) values for stiffness, dissipation and mass lead to an overdamped system (stiffness=" << k << " dissipation=" << disp << " mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return 2.0 * constants::pi / std::sqrt(4.0 * k / mass - mathsFunc::square(disp / mass));
}

Mdouble helpers::computeRestitutionCoefficientFromKAndDispAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass)
{
    if (k <= 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromKandDispAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass) stiffness is not set or has an unexpected value, (stiffness=" << k << ")" << std::endl;
        exit(-1);
    }
    if (disp < 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromKandDispAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass) dissipation is not set or has an unexpected value, (dissipation=" << disp << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromKandDispAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    if (4.0 * mass * k - mathsFunc::square(disp) < 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromKandDispAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass) values for stiffness, dissipation and mass lead to an overdamped system (stiffness=" << k << " dissipation=" << disp << " mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return std::exp(-disp * constants::pi / std::sqrt(4.0 * mass * k - mathsFunc::square(disp)));
}

Mdouble helpers::computeDispFromKAndRestitutionCoefficientAndEffectiveMass(Mdouble k, Mdouble r, Mdouble mass)
{
    if (k <= 0)
    {
        std::cerr << "Error in helpers::computeDispFromKAndRestitutionCoefficientAndEffectiveMass(Mdouble k, Mdouble r, Mdouble mass) stiffness is not set or has an unexpected value, (stiffness=" << k << ")" << std::endl;
        exit(-1);
    }
    if (r < 0)
    {
        std::cerr << "Error in helpers::computeDispFromKAndRestitutionCoefficientAndEffectiveMass(Mdouble k, Mdouble r, Mdouble mass) restitution coefficient is not set or has an unexpected value, (restitution coefficient=" << r << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeDispFromKAndRestitutionCoefficientAndEffectiveMass(Mdouble k, Mdouble r, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return -2.0 * sqrt(mass * k / (constants::sqr_pi + mathsFunc::square(std::log(r)))) * std::log(r);
}

Mdouble helpers::computeCollisionTimeFromKAndRestitutionCoefficientAndEffectiveMass(Mdouble k, Mdouble r, Mdouble mass)
{
    if (k <= 0)
    {
        std::cerr << "Error in helpers::computeCollisionTimeFromKAndRestitutionCoefficientAndEffectiveMass(Mdouble k, Mdouble r, Mdouble mass) stiffness is not set or has an unexpected value, (stiffness=" << k << ")" << std::endl;
        exit(-1);
    }
    if (r < 0)
    {
        std::cerr << "Error in helpers::computeCollisionTimeFromKAndRestitutionCoefficientAndEffectiveMass(Mdouble k, Mdouble r, Mdouble mass) restitution coefficient is not set or has an unexpected value, (restitution coefficient=" << r << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeCollisionTimeFromKAndRestitutionCoefficientAndEffectiveMass(Mdouble k, Mdouble r, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return sqrt(mass / k * (constants::sqr_pi + mathsFunc::square(std::log(r))));
}

Mdouble helpers::computeDispFromKAndCollisionTimeAndEffectiveMass(Mdouble k, Mdouble tc, Mdouble mass)
{
    if (k <= 0)
    {
        std::cerr << "Error in helpers::computeDispFromKAndCollisionTimeAndEffectiveMass(Mdouble k, Mdouble tc, Mdouble mass) stiffness is not set or has an unexpected value, (stiffness=" << k << ")" << std::endl;
        exit(-1);
    }
    if (tc <= 0)
    {
        std::cerr << "Error in helpers::computeDispFromKAndCollisionTimeAndEffectiveMass(Mdouble k, Mdouble tc, Mdouble mass) collision time is not set or has an unexpected value, (collision time=" << tc << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeDispFromKAndCollisionTimeAndEffectiveMass(Mdouble k, Mdouble tc, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    if (mass * k - constants::sqr_pi * mathsFunc::square(mass / tc) < 0)
    {
        std::cerr << "Error in helpers::computeDispFromKAndCollisionTimeAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass) values for stiffness, collision time and mass lead to an overdamped system (stiffness=" << k << " collision time=" << tc << " mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return 2.0 * std::sqrt(mass * k - constants::sqr_pi * mathsFunc::square(mass / tc));
}

Mdouble helpers::computeRestitutionCoefficientFromKAndCollisionTimeAndEffectiveMass(Mdouble k, Mdouble tc, Mdouble mass)
{
    if (k <= 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromKAndCollisionTimeAndEffectiveMass(Mdouble k, Mdouble tc, Mdouble mass) stiffness is not set or has an unexpected value, (stiffness=" << k << ")" << std::endl;
        exit(-1);
    }
    if (tc <= 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromKAndCollisionTimeAndEffectiveMass(Mdouble k, Mdouble tc, Mdouble mass) collision time is not set or has an unexpected value, (collision time=" << tc << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromKAndCollisionTimeAndEffectiveMass(Mdouble k, Mdouble tc, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    if (k / mass * mathsFunc::square(tc) - constants::sqr_pi < 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromKAndCollisionTimeAndEffectiveMass(Mdouble k, Mdouble disp, Mdouble mass) values for stiffness, collision time and mass lead to an overdamped system (stiffness=" << k << " collision time=" << tc << " mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return std::exp(-std::sqrt(k / mass * mathsFunc::square(tc) - constants::sqr_pi));
}

Mdouble helpers::computeDispFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble mass)
{
    if (tc <= 0)
    {
        std::cerr << "Error in helpers::computeDispFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble mass) collision time is not set or has an unexpected value, (collision time=" << tc << ")" << std::endl;
        exit(-1);
    }
    if (r < 0)
    {
        std::cerr << "Error in helpers::computeDispFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble mass) restitution coefficient not set or has an unexpected value, (restitution coefficient=" << r << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeDispFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return -2.0 * mass * std::log(r) / tc;
}

Mdouble helpers::computeKFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble mass)
{
    if (tc <= 0)
    {
        std::cerr << "Error in helpers::computeKFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble mass) collision time is not set or has an unexpected value, (collision time=" << tc << ")" << std::endl;
        exit(-1);
    }
    if (r < 0)
    {
        std::cerr << "Error in helpers::computeKFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble mass) restitution coefficient not set or has an unexpected value, (restitution coefficient=" << r << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeKFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return mass * (mathsFunc::square(constants::pi / tc) + mathsFunc::square(-std::log(r) / tc));
}

Mdouble helpers::computeKFromCollisionTimeAndDispAndEffectiveMass(Mdouble tc, Mdouble disp, Mdouble mass)
{
    if (tc <= 0)
    {
        std::cerr << "Error in helpers::computeKFromCollisionTimeAndDispAndEffectiveMass(Mdouble tc, Mdouble disp, Mdouble mass) collision time is not set or has an unexpected value, (collision time=" << tc << ")" << std::endl;
        exit(-1);
    }
    if (disp < 0)
    {
        std::cerr << "Error in helpers::computeKFromCollisionTimeAndDispAndEffectiveMass(Mdouble tc, Mdouble disp, Mdouble mass) dissipation not set or has an unexpected value, (dissipation=" << disp << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeKFromCollisionTimeAndDispAndEffectiveMass(Mdouble tc, Mdouble disp, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return 0.25 * mathsFunc::square(disp) / mass + constants::sqr_pi * mass / mathsFunc::square(tc);
}

Mdouble helpers::computeRestitutionCoefficientFromCollisionTimeAndDispAndEffectiveMass(Mdouble tc, Mdouble disp, Mdouble mass)
{
    if (tc <= 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromCollisionTimeAndDispAndEffectiveMass(Mdouble tc, Mdouble disp, Mdouble mass) collision time is not set or has an unexpected value, (collision time=" << tc << ")" << std::endl;
        exit(-1);
    }
    if (disp < 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromCollisionTimeAndDispAndEffectiveMass(Mdouble tc, Mdouble disp, Mdouble mass) dissipation not set or has an unexpected value, (dissipation=" << disp << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeRestitutionCoefficientFromCollisionTimeAndDispAndEffectiveMass(Mdouble tc, Mdouble disp, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return std::exp(-0.5 * disp * tc / mass);
}

Mdouble helpers::computeKFromDispAndRestitutionCoefficientAndEffectiveMass(Mdouble disp, Mdouble r, Mdouble mass)
{
    if (disp <= 0)
    {
        std::cerr << "Error in helpers::computeKFromDispAndRestitutionCoefficientAndEffectiveMass(Mdouble disp, Mdouble r, Mdouble mass) dissipation is not set or has an unexpected value, (dissipation time=" << disp << ")" << std::endl;
        exit(-1);
    }
    if (r < 0)
    {
        std::cerr << "Error in helpers::computeKFromDispAndRestitutionCoefficientAndEffectiveMass(Mdouble disp, Mdouble r, Mdouble mass) restitution coefficient not set or has an unexpected value, (restitution coefficient=" << r << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeKFromDispAndRestitutionCoefficientAndEffectiveMass(Mdouble disp, Mdouble r, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return 0.25 * mathsFunc::square(disp)*(constants::sqr_pi / (mathsFunc::square(std::log(r))) + 1.0) / mass;
}

Mdouble helpers::computeCollisionTimeFromDispAndRestitutionCoefficientAndEffectiveMass(Mdouble disp, Mdouble r, Mdouble mass)
{
    if (disp <= 0)
    {
        std::cerr << "Error in helpers::computeCollisionTimeFromDispAndRestitutionCoefficientAndEffectiveMass(Mdouble disp, Mdouble r, Mdouble mass) dissipation is not set or has an unexpected value, (dissipation time=" << disp << ")" << std::endl;
        exit(-1);
    }
    if (r < 0)
    {
        std::cerr << "Error in helpers::computeCollisionTimeFromDispAndRestitutionCoefficientAndEffectiveMass(Mdouble disp, Mdouble r, Mdouble mass) restitution coefficient not set or has an unexpected value, (restitution coefficient=" << r << ")" << std::endl;
        exit(-1);
    }
    if (mass <= 0)
    {
        std::cerr << "Error in helpers::computeCollisionTimeFromDispAndRestitutionCoefficientAndEffectiveMass(Mdouble disp, Mdouble r, Mdouble mass) mass is not set or has an unexpected value (mass=" << mass << ")" << std::endl;
        exit(-1);
    }
    return -2.0 * mass * std::log(r) / disp;
}

helpers::KAndDispAndKtAndDispt helpers::computeDisptFromCollisionTimeAndRestitutionCoefficientAndTangentialRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble beta, Mdouble mass)
{
    helpers::KAndDispAndKtAndDispt ans;
    ans.disp = -2.0 * mass * log(r) / tc;
    ans.k = mass * (mathsFunc::square(constants::pi / tc) + mathsFunc::square(ans.disp / (2.0 * mass)));
    ans.kt = 2.0 / 7.0 * ans.k * (mathsFunc::square(constants::pi) + mathsFunc::square(log(beta))) / (mathsFunc::square(constants::pi) + mathsFunc::square(log(r)));
    if (beta != 0.0)
        ans.dispt = -2 * log(beta) * sqrt(1.0 / 7.0 * mass * ans.kt / (mathsFunc::square(constants::pi) + mathsFunc::square(log(beta))));
    else
        ans.dispt = 2. * sqrt(1.0 / 7.0 * mass * ans.kt);
    return ans;
}

Mdouble helpers::getMaximumVelocity(Mdouble k, Mdouble disp, Mdouble radius, Mdouble mass)
{
    // note: underestimate based on energy argument,
    // Ekin = 2*1/2*m*(v/2)^2 = 1/2*k*(2*r)^2, gives
    // return radius * sqrt(8.0*k/mass);

    // with dissipation, see S. Luding, Collisions & Contacts between two particles, eq 34
    Mdouble w = sqrt(k / mass - mathsFunc::square(disp / (2.0 * mass)));
    Mdouble w0 = sqrt(k / mass);
    Mdouble DispCorrection = exp(-disp / mass / w) * asin(w / w0);
    //std::cout << "DispCorrection" << DispCorrection << std::endl;
    return radius * sqrt(8.0 * k / mass) / DispCorrection;
}

unsigned int helpers::getSaveCountFromNumberOfSavesAndTimeMaxAndTimestep(unsigned int numberOfSaves, Mdouble timeMax, Mdouble timestep)
{
    if (numberOfSaves > 0 && timeMax > 0 && timestep > 0)
    {
        return static_cast<unsigned int>(ceil((timeMax + timestep) / timestep / static_cast<double>(numberOfSaves - 1)));
    }
    else
    {
        logger(ERROR, "[Helpers::getSaveCountFromNumberOfSavesAndTimeMaxAndTimestep()] numberOfSaves: %, timeMax: %, timestep: %", numberOfSaves, timeMax, timestep);
        logger(ERROR, " Arguments need to be positive");
        exit(-1);
    }
}

//seems to be unused, consider taking out \author weinhartt
//unsigned int helpers::getSaveCountFromNumberOfSavesPerTimeUnitAndTimestep(unsigned int numberOfSaves, Mdouble timestep)
//{
//    if (numberOfSaves > 1 && timestep > 0)
//    {
//        return static_cast<unsigned int>(ceil(1.0 / timestep / static_cast<double>(numberOfSaves - 1)));
//    }
//    else
//    {
//        std::cerr << "Error in getSaveCountFromNumberOfSavesPerTimeUnitAndTimestep (" << numberOfSaves << "," << timestep << ")" << std::endl;
//        std::cerr << "Arguments need to be positive" << std::endl;
//        exit(-1);
//    }
//}

void helpers::getLineFromStringStream(std::istream& in, std::stringstream& out)
{
    std::string line_string;
    getline(in, line_string);
    out.str(line_string);
    out.clear();
}

bool helpers::writeToFile(std::string filename, std::string filecontent)
{
    std::fstream file;
    file.open(filename.c_str(), std::ios::out);
    if (file.fail())
    {
        std::cerr << "Error in writeToFile: file could not be opened" << std::endl;
        return false;
    }
    file << filecontent;
    file.close();
    return true;
}

bool helpers::addToFile(std::string filename, std::string filecontent)
{
    std::fstream file;
    file.open(filename.c_str(), std::ios::app);
    if (file.fail())
    {
        std::cerr << "Error in writeToFile: file could not be opened" << std::endl;
        return false;
    }
    file << filecontent;
    file.close();
    return true;
}

bool helpers::fileExists(std::string strFilename)
{
    struct stat stFileInfo;
    bool blnReturn;
    int intStat;

    // Attempt to get the file attributes

    intStat = stat(strFilename.c_str(), &stFileInfo);
    if (intStat == 0)
    {
        // We were able to get the file attributes
        // so the file obviously exists.
        blnReturn = true;
    }
    else
    {
        // We were not able to get the file attributes.
        // This may mean that we don't have permission to
        // access the folder which contains this file. If you
        // need to do that level of checking, lookup the
        // return values of stat which will give you
        // more details on why stat failed.
        blnReturn = false;
    }

    return blnReturn;
}

bool helpers::openFile(std::fstream& file, std::string filename, std::fstream::openmode mode)
{
    file.open(filename.c_str(), mode);
    if (file.fail())
        return false;
    else
        return true;
}

Mdouble helpers::getEffectiveMass(Mdouble mass0, Mdouble mass1)
{
    return mass0 * mass1 / (mass0 + mass1);
}

std::vector<double> helpers::readArrayFromFile(std::string filename, int& n, int& m)
{
    std::fstream file;
    file.open(filename.c_str(), std::ios::in);
    if (file.fail())
    {
        std::cerr << "Error in readArrayFromFile: file could not be opened" << std::endl;
        exit(-1);
    }
    file >> n >> m;
    std::vector<double> v;
    Mdouble val;
    for (int i=0; i<n*m; i++)
    {
        file >> val;
        v.push_back(val);
    }
    file.close();
    return v;
}

void helpers::loadingTest(const ParticleSpecies* species, Mdouble displacement, Mdouble velocity, Mdouble radius, std::string name)
{
    class LoadingTest : public DPMBase {
        const ParticleSpecies* species;
        Mdouble displacement;
        Mdouble velocity;
        Mdouble radius;
    public:
        //public variables
        LoadingTest (const ParticleSpecies* species, Mdouble displacement, Mdouble velocity, Mdouble radius)
        : species(species), displacement(displacement), velocity(velocity), radius(radius) {}

        void setupInitialConditions() override
        {
            //setName("LoadingTest"+species->getName());
            setTimeMax(2.0*displacement/velocity);
            setTimeStep(2e-3*getTimeMax());
            setSaveCount(1);
            setFileType(FileType::NO_FILE);
            fStatFile.setFileType(FileType::ONE_FILE);

            setMax({radius,radius,radius+radius});
            setMin({-radius,-radius,0});
            setSystemDimensions(3);
            setParticleDimensions(3);

            speciesHandler.copyAndAddObject(*species);

            BaseParticle p;
            p.setSpecies(speciesHandler.getObject(0));
            p.setRadius(radius);
            p.setPosition({0,0,radius});
            particleHandler.copyAndAddObject(p);

            InfiniteWall w;
            w.setSpecies(speciesHandler.getObject(0));
            w.set(Vec3D(0,0,-1),Vec3D(0.0,0.0,0.0));
            wallHandler.copyAndAddObject(w);
        }

        void actionsBeforeTimeStep() override {
            BaseParticle* p = particleHandler.getLastObject();
            assert(p);
            p->setAngularVelocity({0,0,0});

            //Moving particle normally into surface
            if (getTime() <= displacement/velocity) {
                p->setVelocity({0,0,velocity});
                p->setPosition({0,0,radius-velocity*getTime()});
            } else {
                p->setVelocity({0,0,-velocity});
                p->setPosition({0,0,radius-displacement-displacement+velocity*getTime()});
            }
        }
    } test(species, displacement, velocity, radius);
    test.setName(name);
    test.solve();
    writeToFile(test.getName()+".gnu","plot '"+test.getName()+".fstat' u 7:9 w lp");
    logger(INFO,"finished loading test: run 'gnuplot %.gnu' to view output",test.getName());
}

void helpers::normalAndTangentialLoadingTest(const ParticleSpecies *species, Mdouble displacement,
                                             Mdouble tangentialDisplacement, Mdouble velocity, Mdouble radius,
                                             std::string name)
{
    class LoadingTest : public DPMBase {
        const ParticleSpecies* species;
        Mdouble displacement;
        Mdouble tangentialDisplacement;
        Mdouble velocity;
        Mdouble radius;
    public:
        //public variables
        LoadingTest (const ParticleSpecies* species, Mdouble displacement, Mdouble tangentialDisplacement, Mdouble velocity, Mdouble radius)
         : species(species), displacement(displacement), tangentialDisplacement(tangentialDisplacement), velocity(velocity), radius(radius) {}

        void setupInitialConditions() override
        {
            //setName("TangentialLoadingTest"+species->getName());
            setTimeMax(4.0*tangentialDisplacement/velocity);
            setTimeStep(2e-3*getTimeMax());
            setSaveCount(1);
            setFileType(FileType::NO_FILE);
            fStatFile.setFileType(FileType::ONE_FILE);

            setMax({radius,radius,radius+radius});
            setMin({-radius,-radius,0});
            setSystemDimensions(3);
            setParticleDimensions(3);

            speciesHandler.copyAndAddObject(*species);

            BaseParticle p;
            p.setSpecies(speciesHandler.getObject(0));
            p.setRadius(radius);
            p.setPosition({0,0,radius-displacement});
            particleHandler.copyAndAddObject(p);

            InfiniteWall w;
            w.setSpecies(speciesHandler.getObject(0));
            w.set(Vec3D(0,0,-1),Vec3D(0.0,0.0,0.0));
            wallHandler.copyAndAddObject(w);
        }

        void actionsBeforeTimeStep() override {
            BaseParticle* p = particleHandler.getLastObject();
            assert(p);
            p->setAngularVelocity({0,0,0});

            //Moving particle normally into surface
            Mdouble time = getTime() / (tangentialDisplacement/velocity);
            if (time <= 1.0) {
                p->setVelocity({velocity,0,0});
                p->setPosition({velocity*getTime(),0,radius-displacement});
            } else if (time <= 3.0) {
                p->setVelocity({-velocity,0,0});
                p->setPosition({2.0*tangentialDisplacement-velocity*getTime(),0,radius-displacement});
            } else {
                p->setVelocity({velocity,0,0});
                p->setPosition({-4.0*tangentialDisplacement+velocity*getTime(),0,radius-displacement});
            }
        }

    } test(species, displacement, tangentialDisplacement, velocity, radius);
    test.setName(name);
    test.solve();
    writeToFile(test.getName()+".gnu","plot '"+test.getName()+".fstat' u 8:($10*$14) w lp");
    logger(INFO,"finished tangential loading test: run 'gnuplot %.gnu' to view output",test.getName());
}

void helpers::objectivenessTest(const ParticleSpecies* species, Mdouble displacement, Mdouble tangentialDisplacement, Mdouble velocity, Mdouble radius, std::string name)
{
    class ObjectivenessTest : public DPMBase {
        const ParticleSpecies* species;
        Mdouble displacement;
        Mdouble tangentialDisplacement;
        Mdouble velocity;
        Mdouble radius;
    public:
        //public variables
        ObjectivenessTest (const ParticleSpecies* species, Mdouble displacement, Mdouble tangentialDisplacement, Mdouble velocity, Mdouble radius)
         : species(species), displacement(displacement), tangentialDisplacement(tangentialDisplacement), velocity(velocity), radius(radius) {}

        void setupInitialConditions() override
        {
            //setName("ObjectivenessTest"+species->getName());
            setTimeMax((tangentialDisplacement+0.5*constants::pi*radius)/velocity);
            setTimeStep(1e-4*getTimeMax());
            setSaveCount(20);
            setFileType(FileType::NO_FILE);
            dataFile.setFileType(FileType::ONE_FILE);
            fStatFile.setFileType(FileType::ONE_FILE);

            setMax(radius*Vec3D(2,2,2));
            setMin(radius*Vec3D(-2,-2,-2));
            setSystemDimensions(2);
            setParticleDimensions(3);

            speciesHandler.copyAndAddObject(*species);

            BaseParticle p;
            p.setSpecies(speciesHandler.getObject(0));
            p.setRadius(radius);
            p.setPosition({0,radius-displacement,0});
            particleHandler.copyAndAddObject(p);
            p.setPosition({0,-radius+displacement,0});
            particleHandler.copyAndAddObject(p);
        }

        void actionsBeforeTimeStep() override {
            BaseParticle* p = particleHandler.getObject(0);
            BaseParticle* q = particleHandler.getLastObject();
            assert(p);
            assert(q);

            //Moving particle normally into surface
            if (getTime() <= tangentialDisplacement/velocity) {
                p->setAngularVelocity({0,0,0});
                p->setVelocity({velocity,0,0});
                p->setPosition({-tangentialDisplacement+velocity*getTime(),radius-displacement,0});
                q->setAngularVelocity({0,0,0});
                q->setVelocity({-velocity,0,0});
                q->setPosition({tangentialDisplacement-velocity*getTime(),-radius+displacement,0});
            } else {
                Mdouble angle = velocity/(radius-displacement) * (getTime()-tangentialDisplacement/velocity);
                Mdouble s = sin(angle);
                Mdouble c = cos(angle);
                p->setAngularVelocity(velocity/(radius-displacement)*Vec3D(0,0,-1));
                //p->setAngularVelocity(Vec3D(0,0,0));
                //p->setOrientation({1,0,0,0});
                p->setVelocity(velocity*Vec3D(c,-s,0));
                //p->setVelocity(Vec3D(0,0,0));
                p->setPosition((radius-displacement)*Vec3D(s,c,0));
                q->setAngularVelocity(-p->getAngularVelocity());
                q->setOrientation(-p->getOrientation());
                q->setVelocity(-p->getVelocity());
                q->setPosition(-p->getPosition());
            }
        }

    } test (species, displacement, tangentialDisplacement, velocity, radius);
    test.setName(name);
    test.solve();
    writeToFile(test.getName()+".gnu","set size ratio -1; plot '"+test.getName()+".fstat' u 14:15 every 2 w lp");
    logger(INFO,"finished objectiveness test: run 'gnuplot %.gnu' to view output",test.getName());
}

bool helpers::compare(std::istream& is,std::string s)
{
    // Get current position
    //check if the next line starts with 'interactionFile'; otherwise, skip interaction
    int len = is.tellg();
    std::string dummy;
    is >> dummy;
    if (dummy.compare(s) != 0)
    {
        is.seekg(len, std::ios_base::beg);
        return false;
        logger(INFO, "helpers::compare: Next stream value (%) is not %", dummy,s);
    }
    return true;
}

void helpers::check(double real, double ideal, double error, std::string errorMessage) {
    logger.assert_always(mathsFunc::isEqual(real, ideal, error),
                         errorMessage + ": % (should be %) ",real, ideal);
}