Codes for tutorials

Particle motion in outer space (code)

Goto tutorial T1: Particle motion in outer space

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T1
*/
 
#include <Mercury3D.h>
#include <Species/Species.h>
#include <Species/LinearViscoelasticSpecies.h>
 
class Tutorial1 : public Mercury3D 
{
 
public:
 
    void setupInitialConditions() override {
        // create a new particle
        SphericalParticle p0;
        // set species (material type)
        p0.setSpecies(speciesHandler.getObject(0));
        // set particle radius, position, velocity
        p0.setRadius(0.01);
        p0.setPosition(Vec3D(0.1 , 0.1, 0.1 ));
        p0.setVelocity(Vec3D(0.1, 0.0, 0.0));
        // pass particle to the particle handler (which contains all particles in the simulation)
        particleHandler.copyAndAddObject(p0);
    }
};
 
int main(int argc, char* argv[])
{
    // Instantiate an object of class Tutorial1
    Tutorial1 problem;
 
    // Problem setup: set file name, gravity, domain size, simulation time
    problem.setName("Tutorial1");
    problem.setGravity(Vec3D(0.0, 0.0, 0.0));
    problem.setXMax(1.0);
    problem.setYMax(1.0);
    problem.setZMax(1.0);
    problem.setTimeMax(1.0);
 
    // Set the species (material properties such as density and stiffness) of particles and walls
    // (in this case, both are assumed to consist of the same material)
    LinearViscoelasticSpecies species;
    species.setDensity(2500.0); //sets the species type_0 density in kg/m3
    // The normal spring stiffness and normal dissipation is computed such that
    // collision time tc=0.005 and restitution coefficient rc=1.0
    species.setStiffness(258.5);//sets the spring stiffness in kN/m.
    species.setDissipation(0.0); //sets the dissipation.
    problem.speciesHandler.copyAndAddObject(species);
 
    // Define what output gets written
    // number of time steps skipped between saves (i.e. every 10-th time step is written to file)
    problem.setSaveCount(50);
    // creates file with particle positions, velocities, radii for multiple time steps (for plotting)
    problem.dataFile.setFileType(FileType::ONE_FILE);
    // file with contact forces, overlaps for multiple time steps (for plotting)
    problem.fStatFile.setFileType(FileType::NO_FILE);
    // file with all parameters of the last saved time step (for restarting)
    problem.restartFile.setFileType(FileType::ONE_FILE);
    // writes global properties like kinetic/potential energy and center of mass (for quick analysis)
    problem.eneFile.setFileType(FileType::NO_FILE);
 
    // whether the wall geometry is written into a vtu file (either once initially or for several time steps)
    problem.wallHandler.setWriteVTK(FileType::ONE_FILE);
    // whether the particle positions are written into a vtu file
    problem.particleHandler.setWriteVTK(true);
 
    // set a time step 1/50th of collision time
    problem.setTimeStep(0.005 / 50.0);
    // start the solver (calls setupInitialConditions and initiates time loop)
    problem.solve(argc, argv);
 
    return 0;
}

Return to tutorial T1: Particle motion in outer space

Particle motion on earth (code)

Goto tutorial T2: Particle motion on earth

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Tutorial 2
 
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T2
*/
 
#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
 
class Tutorial2 : public Mercury3D
{
public:
    
    void setupInitialConditions() override {
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setRadius(0.05);
        p0.setPosition(Vec3D(0.5 * getXMax(), 0.5 * getYMax(), getZMax()));
        p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
        particleHandler.copyAndAddObject(p0);
    }
    
};
 
int main(int argc, char* argv[])
{
    // Problem setup
    Tutorial2 problem;
 
    problem.setName("Tutorial2");
    problem.setGravity(Vec3D(0.0, 0.0, -9.81));
    problem.setXMax(1.0);
    problem.setYMax(1.0);
    problem.setZMax(5.0);
    problem.setTimeMax(1.5);
 
    // The normal spring stiffness and normal dissipation is computed and set as
    // For collision time tc=0.005 and restitution coefficient rc=1.0,
    LinearViscoelasticSpecies species;
    species.setDensity(2500.0); //sets the species type_0 density
    species.setStiffness(258.5);//sets the spring stiffness.
    species.setDissipation(0.0); //sets the dissipation.
    problem.speciesHandler.copyAndAddObject(species);
 
    problem.setSaveCount(50);
    problem.dataFile.setFileType(FileType::ONE_FILE);
    problem.restartFile.setFileType(FileType::ONE_FILE);
    problem.fStatFile.setFileType(FileType::NO_FILE);
    problem.eneFile.setFileType(FileType::NO_FILE);
 
    // whether the wall geometry is written into a vtu file (either once initially or for several time steps)
    problem.wallHandler.setWriteVTK(FileType::ONE_FILE);
    // whether the particle positions are written into a vtu file
    problem.particleHandler.setWriteVTK(true);
 
 
    problem.setTimeStep(.005 / 50.0);// (collision time)/50.0
    problem.solve(argc, argv);
    
    return 0;
}

Return to tutorial T2: Particle motion on earth

Bouncing ball - elastic (code)

Goto tutorial T3: Bouncing ball (elastic)

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Tutorial 3
 
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T3
*/
 
#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
#include <Walls/InfiniteWall.h>
 
class Tutorial3 : public Mercury3D
{
public:
 
    void setupInitialConditions() override {
        // add a particle of 1cm diameter at height zMax
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setRadius(0.005);
        p0.setPosition(Vec3D(0.5 * getXMax(), 0.5 * getYMax(), getZMax()));
        p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
        particleHandler.copyAndAddObject(p0);
 
        // add a bottom wall at zMin
        InfiniteWall w0;
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(0.0, 0.0, -1.0), Vec3D(0, 0, getZMin()));
        wallHandler.copyAndAddObject(w0);
    }
    
};
 
int main(int argc, char* argv[])
{
    
    // Problem setup
    Tutorial3 problem;
 
    // name the output files
    problem.setName("Tutorial3");
    // remove old output files with the same name
    problem.removeOldFiles();
    // set gravivational acceleration
    problem.setGravity(Vec3D(0.0, 0.0, -9.81));
    // set domain size (xMin = yMin = zMin = 0 by default)
    problem.setXMax(0.5);
    problem.setYMax(0.5);
    problem.setZMax(0.5);
    problem.setTimeMax(2.0);
 
    // Sets a linear spring-damper contact law.
    // The normal spring stiffness and normal dissipation is computed and set as
    // For collision time tc=0.005 and restitution coefficient rc=1.0,
    auto species = problem.speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
    species->setDensity(2500.0);
    double mass = species->getMassFromRadius(0.005);
    double collisionTime = 0.005;
    double restitution = 1.0;
    species->setCollisionTimeAndRestitutionCoefficient(collisionTime,restitution,mass);
    logger(INFO, "Stiffness %, dissipation %", species->getStiffness(), species->getDissipation());
 
    problem.setSaveCount(25);
 
    // Output vtk files for ParaView visualisation
    problem.wallHandler.setWriteVTK(FileType::ONE_FILE);
    problem.particleHandler.setWriteVTK(true);
 
    // Sets time step to 1/50th of the collision time
    problem.setTimeStep(0.005 / 50.0);
    // Start the solver (calls setupInitialConditions and initiates time loop)
    problem.solve(argc, argv);
    return 0;
}

Return to tutorial T3: Bouncing ball (elastic)

Bouncing ball - inelastic (code)

Goto tutorial T4: Bouncing ball with dissipation (inelastic)

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Tutorial 4
 
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T4
*/
 
#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
#include <Walls/InfiniteWall.h>
 
class Tutorial4 : public Mercury3D
{
public:
    
    void setupInitialConditions() override {
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setRadius(0.005);
        p0.setPosition(Vec3D(0.5 * getXMax(), 0.5 * getYMax(), getZMax()));
        p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
        particleHandler.copyAndAddObject(p0);
        
        InfiniteWall w0;
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(0.0, 0.0, -1.0), Vec3D(0.0, 0.0, getZMin()));
        wallHandler.copyAndAddObject(w0);
    }
 
};
 
int main(int argc, char* argv[])
{
    
    // Problem setup
    Tutorial4 problem;
    
    problem.setName("Tutorial4");
    problem.setGravity(Vec3D(0.0, 0.0, -9.81));
    problem.setXMax(0.25);
    problem.setYMax(0.25);
    problem.setZMax(0.25);
    problem.setTimeMax(5.0);
 
    // The normal spring stiffness and normal dissipation is computed and set as
    // For collision time tc=0.005 and restitution coefficient rc=0.88,
    LinearViscoelasticSpecies species;
    species.setDensity(2500.0); //sets the species type_0 density
    species.setStiffness(258.5);//sets the spring stiffness.
    species.setDissipation(0.0334); //sets the dissipation.
    problem.speciesHandler.copyAndAddObject(species);
    
    problem.setSaveCount(10);
    problem.setFileType(FileType::ONE_FILE);
    problem.fStatFile.setFileType(FileType::NO_FILE);
 
    // TODO: Write vtk files (exercise 5)
 
    //time integration parameters
    problem.setTimeStep(0.005 / 50.0); // (collision time)/50.0
    problem.solve(argc, argv);
    return 0;
}

Return to tutorial T4: Bouncing ball with dissipation (inelastic)

Elastic collision - 2 particles (code)

Goto tutorial T5: Elastic collision (2 particles)

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Tutorial 5
 
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T5
*/
 
#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
 
class Tutorial5 : public Mercury3D
{
public:
    
    void setupInitialConditions() override {
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        
        p0.setRadius(0.005); //particle-1 radius
        p0.setPosition(Vec3D(0.25 * getXMax(), 0.5 * getYMax(), 0.5 * getZMax()));
        p0.setVelocity(Vec3D(0.25, 0.0, 0.0));
        particleHandler.copyAndAddObject(p0); // 1st particle created
        
        p0.setRadius(0.005); // particle-2 radius
        p0.setPosition(Vec3D(0.75 * getXMax(), 0.5 * getYMax(), 0.5 * getZMax()));
        p0.setVelocity(Vec3D(-0.25, 0.0, 0.0));
        particleHandler.copyAndAddObject(p0); // 2nd particle created
    }
    
};
 
int main(int argc, char* argv[])
{
    
    // Problem setup
    Tutorial5 problem;
    
    problem.setName("Tutorial5");
    problem.setGravity(Vec3D(0.0, 0.0, 0.0));
    problem.setXMax(0.5);
    problem.setYMax(0.25);
    problem.setZMax(0.5);
    problem.setTimeMax(2.0);
 
    // The normal spring stiffness and normal dissipation is computed and set as
    // For collision time tc=0.005 and restitution coefficient rc=1.0,
    LinearViscoelasticSpecies species;
    species.setDensity(2500.0); //sets the species type_0 density
    species.setStiffness(258.5);//sets the spring stiffness.
    species.setDissipation(0.0); //sets the dissipation.
    problem.speciesHandler.copyAndAddObject(species);
    
    problem.setSaveCount(10);
    problem.dataFile.setFileType(FileType::ONE_FILE);
    problem.restartFile.setFileType(FileType::ONE_FILE);
    problem.fStatFile.setFileType(FileType::NO_FILE);
    problem.eneFile.setFileType(FileType::NO_FILE);
 
    // TODO: Write vtk files (exercise 1)
 
    problem.setTimeStep(.005 / 50.0); // (collision time)/50.0
    problem.solve(argc, argv);
    
    return 0;
}

Return to tutorial T5: Elastic collision (2 particles)

Elastic collisions with periodic boundaries (code)

Goto tutorial T6: Elastic collisions with periodic boundaries

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Tutorial 6
 
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T6
*/
 
#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
#include <Boundaries/PeriodicBoundary.h>
 
class Tutorial6 : public Mercury3D
{
public:
    
    void setupInitialConditions() override {
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setRadius(0.015);//particle-1
        p0.setPosition(Vec3D(0.15 * getXMax(), 0.5 * getYMax(), 0.5 * getZMax()));
        p0.setVelocity(Vec3D(0.25, 0.0, 0.0));
        particleHandler.copyAndAddObject(p0);
        
        p0.setRadius(0.005);// particle-2
        p0.setPosition(Vec3D(0.75 * getXMax(), 0.5 * getYMax(), 0.5 * getZMax()));
        p0.setVelocity(Vec3D(-0.25, 0.0, 0.0));
        particleHandler.copyAndAddObject(p0);
 
        PeriodicBoundary b0;
        b0.set(Vec3D(1.0, 0.0, 0.0), getXMin(), getXMax());
        boundaryHandler.copyAndAddObject(b0);
    }
    
};
 
int main(int argc, char* argv[])
{
    
    // Problem setup
    Tutorial6 problem; // instantiate an object of class Tutorial 6
    
    problem.setName("Tutorial6");
    problem.setGravity(Vec3D(0.0, 0.0, 0.0));
    problem.setXMax(0.5);
    problem.setYMax(0.25);
    problem.setZMax(0.5);
    problem.setTimeMax(5.0);
 
    // The normal spring stiffness and normal dissipation is computed and set as
    // For collision time tc=0.005 and restitution coefficient rc=1.0,
    LinearViscoelasticSpecies species;
    species.setDensity(2500.0); //sets the species type_0 density
    species.setStiffness(258.5);//sets the spring stiffness.
    species.setDissipation(0.0); //sets the dissipation.
    problem.speciesHandler.copyAndAddObject(species);
 
    
    problem.setSaveCount(10);
    problem.dataFile.setFileType(FileType::ONE_FILE);
    problem.restartFile.setFileType(FileType::ONE_FILE);
    problem.fStatFile.setFileType(FileType::NO_FILE);
    problem.eneFile.setFileType(FileType::NO_FILE);
 
    // Write vtk files
    problem.wallHandler.setWriteVTK(FileType::ONE_FILE);
    problem.particleHandler.setWriteVTK(true);
    
    problem.setTimeStep(0.005 / 50.0);// (collision time)/50.0
    problem.solve(argc, argv);
    
    return 0;
}

Return to tutorial T6: Elastic collisions with periodic boundaries

Motion of a particle in a two dimensional box (code)

Goto tutorial T7: Motion of a particle in a two dimensional (2D) box

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Tutorial 7
 
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T7
*/
 
#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
#include <Walls/InfiniteWall.h>
 
class Tutorial7 : public Mercury3D
{
public:
    
    void setupInitialConditions() override {
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setRadius(0.005);
        p0.setPosition(Vec3D(0.5 * getXMax(), 0.5 * getYMax(), 0.0));
        p0.setVelocity(Vec3D(1.2, 1.3, 0.0));
        particleHandler.copyAndAddObject(p0);
        
        InfiniteWall w0;
        
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(1.0, 0.0, 0.0), Vec3D(getXMax(), 0.0, 0.0));
        wallHandler.copyAndAddObject(w0);
        
        w0.set(Vec3D(-1.0, 0.0, 0.0), Vec3D(getXMin(), 0.0, 0.0));
        wallHandler.copyAndAddObject(w0);
        
        w0.set(Vec3D(0.0, 1.0, 0.0), Vec3D(0.0, getYMax(), 0.0));
        wallHandler.copyAndAddObject(w0);
        
        w0.set(Vec3D(0.0, -1.0, 0.0), Vec3D(0.0, getYMin(), 0.0));
        wallHandler.copyAndAddObject(w0);
    }
    
};
 
int main(int argc, char* argv[])
{
    // Problem setup
    Tutorial7 problem; // instantiate an object of class Tutorial 7
    
    problem.setName("Tutorial7");
    problem.setGravity(Vec3D(0.0, 0.0, 0.0));
    problem.setXMax(1);
    problem.setYMax(0.5);
    problem.setZMax(1.0);
    problem.setTimeMax(2.5);
 
    // The normal spring stiffness and normal dissipation is computed and set as
    // For collision time tc=0.005 and restitution coefficient rc=1.0,
    LinearViscoelasticSpecies species;
    species.setDensity(2500.0); //sets the species type_0 density
    species.setStiffness(258.5);//sets the spring stiffness.
    species.setDissipation(0.0); //sets the dissipation.
    problem.speciesHandler.copyAndAddObject(species);
    
    problem.setSaveCount(50);
    problem.dataFile.setFileType(FileType::ONE_FILE);
    problem.restartFile.setFileType(FileType::ONE_FILE);
    problem.fStatFile.setFileType(FileType::NO_FILE);
    problem.eneFile.setFileType(FileType::NO_FILE);
 
    problem.wallHandler.setWriteVTK(FileType::ONE_FILE);
    problem.particleHandler.setWriteVTK(true);
 
    problem.setTimeStep(0.005 / 50.0);// (collision time)/50.0
    problem.solve(argc, argv);
    
    return 0;
}

Return to tutorial T7: Motion of a particle in a two dimensional (2D) box

Motion of a particle in a box with an obstacle (code)

Goto tutorial T8: Motion of a particle in a box with an obstacle

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Tutorial 8
 
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T8
*/
 
#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
#include <Walls/InfiniteWall.h>
#include <Walls/IntersectionOfWalls.h>
 
class Tutorial8 : public Mercury3D
{
public:
    
    void setupInitialConditions() override {
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setRadius(0.005);
        p0.setPosition(Vec3D(0.15 * getXMax(), 0.3335 * getYMax(), 0.0));
        p0.setVelocity(Vec3D(1.2, 0.2, 0.0));
        particleHandler.copyAndAddObject(p0);
 
        InfiniteWall w0;
        
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(1.0, 0.0, 0.0), Vec3D(getXMax(), 0.0, 0.0));
        wallHandler.copyAndAddObject(w0);
        
        w0.set(Vec3D(-1.0, 0.0, 0.0), Vec3D(getXMin(), 0.0, 0.0));
        wallHandler.copyAndAddObject(w0);
        
        w0.set(Vec3D(0.0, 1.0, 0.0), Vec3D(0.0, getYMax(), 0.0));
        wallHandler.copyAndAddObject(w0);
        
        w0.set(Vec3D(0.0, -1.0, 0.0), Vec3D(0.0, getYMin(), 0.0));
        wallHandler.copyAndAddObject(w0);
 
        IntersectionOfWalls w1;
        w1.setSpecies(speciesHandler.getObject(0));
        
        w1.addObject(Vec3D(-1.0, 0.0, 0.0), Vec3D(0.75 * getXMax(), 0.0, 0.0));
        w1.addObject(Vec3D(1.0, 0.0, 0.0), Vec3D(0.25 * getXMax(), 0.0, 0.0));
        w1.addObject(Vec3D(0.0, -1.0, 0.0), Vec3D(0.0, 0.75 * getYMax(), 0.0));
        w1.addObject(Vec3D(0.0, 1.0, 0.0), Vec3D(0.0, 0.25 * getYMax(), 0.0));
        wallHandler.copyAndAddObject(w1);
    }
    
};
 
int main(int argc, char* argv[])
{
    
    // Problem setup
    Tutorial8 problem; // instantiate an object of class Tutorial 8
    
    problem.setName("Tutorial8");
    problem.setGravity(Vec3D(0.0, 0.0, 0.0));
    problem.setXMax(0.5);
    problem.setYMax(0.5);
    problem.setZMax(0.5);
    problem.setTimeMax(5.0);
 
    // The normal spring stiffness and normal dissipation is computed and set as
    // For collision time tc=0.005 and restitution coefficient rc=1.0,
    LinearViscoelasticSpecies species;
    species.setDensity(2500.0); //sets the species type_0 density
    species.setStiffness(258.5);//sets the spring stiffness.
    species.setDissipation(0.0); //sets the dissipation.
    problem.speciesHandler.copyAndAddObject(species);
 
    
    problem.setSaveCount(50);
    problem.dataFile.setFileType(FileType::ONE_FILE);
    problem.restartFile.setFileType(FileType::ONE_FILE);
    problem.fStatFile.setFileType(FileType::NO_FILE);
    problem.eneFile.setFileType(FileType::NO_FILE);
 
    problem.wallHandler.setWriteVTK(FileType::ONE_FILE);
    problem.particleHandler.setWriteVTK(true);
 
    problem.setTimeStep(.005 / 50.0); // (collision time)/50.0
    problem.solve(argc, argv);
    
    return 0;
}

Return to tutorial T8: Motion of a particle in a box with an obstacle

Motion of a ball over an inclined plane (code)

Goto tutorial T9: Motion of a ball over an inclined plane (Sliding + Rolling)

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Tutorial 9
 
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T9
*/
 
#include <Species/LinearViscoelasticFrictionSpecies.h>
#include <Mercury3D.h>
#include <Walls/InfiniteWall.h>
 
class Tutorial9 : public Mercury3D
{
public:
    
    void setupInitialConditions() override {
        SphericalParticle p0;
 
        //sets the particle to species type-1
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setRadius(0.005);
        p0.setPosition(Vec3D(0.05 * getXMax(), 0.25 * getYMax(), getZMin() + p0.getRadius()));
        p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
        particleHandler.copyAndAddObject(p0);
        
        //sets the particle to species type-2
        p0.setSpecies(speciesHandler.getObject(1));
        p0.setRadius(0.005);
        p0.setPosition(Vec3D(0.05 * getXMax(), 0.5 * getYMax(), getZMin() + p0.getRadius()));
        p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
        particleHandler.copyAndAddObject(p0);
        
        //sets the particle to species type-3
        p0.setSpecies(speciesHandler.getObject(2));
        p0.setRadius(0.005);
        p0.setPosition(Vec3D(0.05 * getXMax(), 0.75 * getYMax(), getZMin() + p0.getRadius()));
        p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
        particleHandler.copyAndAddObject(p0);
 
        InfiniteWall w0;
        
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(0.0, 0.0, -1.0), Vec3D(0.0, 0.0, getZMin()));
        wallHandler.copyAndAddObject(w0);
    }
};
 
int main(int argc, char* argv[])
{
 
    // Instantiate an object of class Tutorial 9
    Tutorial9 problem;
 
    // Problem setup
    problem.setName("Tutorial9");
    problem.removeOldFiles();
    double angle = constants::pi / 180.0 * 20.0;
    problem.setGravity(Vec3D(sin(angle), 0.0, -cos(angle)) * 9.81);
    problem.setXMax(0.3);
    problem.setYMax(0.3);
    problem.setZMax(0.05);
    problem.setTimeMax(0.4);
 
    // Set the species (material properties such as density and stiffness) of particles and walls
    // (in this case, particle-0 and the wall are made from species-0, particle-1 belongs to species-1, particle-2 belongs to species-2)
 
    // Species-0 properties
    LinearViscoelasticFrictionSpecies species0;
    // The normal spring stiffness and normal dissipation is computed and set as
    // For collision time tc=0.005 and restitution coefficient rc=0.88,
    species0.setDensity(2500.0);//sets the species type_0 density
    species0.setStiffness(259.018);//sets the spring stiffness.
    species0.setDissipation(0.0334);//sets the dissipation.
    species0.setSlidingStiffness(2.0 / 7.0 * species0.getStiffness());
    species0.setRollingStiffness(2.0 / 5.0 * species0.getStiffness());
    species0.setSlidingFrictionCoefficient(0.0);
    species0.setRollingFrictionCoefficient(0.0);
    auto ptrToSp0=problem.speciesHandler.copyAndAddObject(species0);
 
    // Species-1 properties
    LinearViscoelasticFrictionSpecies species1;
    species1.setDensity(2500.0);//sets the species type-1 density
    species1.setStiffness(259.018);//sets the spring stiffness
    species1.setDissipation(0.0334);//sets the dissipation
    species1.setSlidingStiffness(2.0 / 7.0 * species1.getStiffness());
    species1.setRollingStiffness(2.0 / 5.0 * species1.getStiffness());
    species1.setSlidingFrictionCoefficient(0.5);
    species1.setRollingFrictionCoefficient(0.0);
    auto ptrToSp1=problem.speciesHandler.copyAndAddObject(species1);
 
    // Properties of contacts between species-0 and species-1
    // (no new species is defined since the mixed species is automatically created)
    auto species01 = problem.speciesHandler.getMixedObject(ptrToSp0,ptrToSp1);
    species01->setStiffness(259.018);//sets the spring stiffness
    species01->setDissipation(0.0334);//sets the dissipation
    species01->setSlidingStiffness(2.0 / 7.0 * species01->getStiffness());
    species01->setRollingStiffness(2.0 / 5.0 * species01->getStiffness());
    species01->setSlidingFrictionCoefficient(0.5);
    species01->setRollingFrictionCoefficient(0.0);
 
    // Species 2 properties
    LinearViscoelasticFrictionSpecies species2;
    species2.setDensity(2500.0);//sets the species type-2 density
    species2.setStiffness(258.5);//sets the spring stiffness
    species2.setDissipation(0.0);//sets the dissipation
    species2.setSlidingStiffness(2.0 / 7.0 * species2.getStiffness());
    species2.setRollingStiffness(2.0 / 5.0 * species2.getStiffness());
    species2.setSlidingFrictionCoefficient(0.5);
    species2.setRollingFrictionCoefficient(0.5);
    auto ptrToSp2 = problem.speciesHandler.copyAndAddObject(species2);
 
    // Properties of contacts between species-0 and species-2
    auto species02 = problem.speciesHandler.getMixedObject(ptrToSp0, ptrToSp2);
    species02->setStiffness(259.018);//sets the stiffness
    species02->setDissipation(0.0334);//sets the dissipation
    species02->setSlidingStiffness(2.0 / 7.0 * species02->getStiffness());
    species02->setRollingStiffness(2.0 / 5.0 * species02->getStiffness());
    species02->setSlidingFrictionCoefficient(0.5);
    species02->setRollingFrictionCoefficient(0.5);
 
    // Note: Properties of contacts between species-0 and species-1 are never defined, since no contacts between those
    // two materials occurs in this simulation; note: if no properties are defined, then default properties are assumed,
    // mostly harmonic means of the properties of both individual species
 
    //Output
    problem.setSaveCount(25);
    problem.dataFile.setFileType(FileType::ONE_FILE);
    problem.restartFile.setFileType(FileType::ONE_FILE);
    problem.fStatFile.setFileType(FileType::ONE_FILE);
    problem.eneFile.setFileType(FileType::NO_FILE);
 
    // whether the wall geometry is written into a vtu file (either once initially or for several time steps)
    problem.wallHandler.setWriteVTK(FileType::ONE_FILE);
    // whether the particle positions are written into a vtu file
    problem.particleHandler.setWriteVTK(true);
 
    // set a time step 1/50th of collision time
    problem.setTimeStep(0.005 / 50.0);
    // start the solver (calls setupInitialConditions and initiates time loop)
    problem.solve(argc, argv);
    return 0;
}

Return to tutorial T9: Motion of a ball over an inclined plane (Sliding + Rolling)

Loading data restart(code)

Goto tutorial T10: How to load restart and data files

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Tutorial 10: This tutorial shows how to load restart and data files
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T10
*/
 
#include "Mercury3D.h"
#include <Math/Helpers.h>
#include "Walls/IntersectionOfWalls.h"
#include "Walls/AxisymmetricIntersectionOfWalls.h"
#include "Species/LinearViscoelasticSpecies.h"
 
 
int main(int argc, char* argv[])
{
 
    //writeToFile is used here to create a restart and a data file, which will be loaded below.
    helpers::writeToFile("Tutorial10.ini.restart",
        "restart_version 1.0 name Tutorial10\n"
        "dataFile name Tutorial10.data fileType ONE_FILE saveCount 10 counter 0 nextSavedTimeStep 0\n"
        "fStatFile name Tutorial10.fstat fileType NO_FILE saveCount 10 counter 0 nextSavedTimeStep 0\n"
        "eneFile name Tutorial10.ene fileType ONE_FILE saveCount 10 counter 0 nextSavedTimeStep 0\n"
        "restartFile name Tutorial10.restart fileType ONE_FILE saveCount 10 counter 0 nextSavedTimeStep 0\n"
        "statFile name Tutorial10.stat fileType ONE_FILE saveCount 10 counter 0 nextSavedTimeStep 0\n"
        "xMin 0 xMax 2 yMin 0 yMax 2 zMin 0 zMax 2\n"
        "timeStep 1e-03 time 0 ntimeSteps 0 timeMax 10\n"
        "systemDimensions 3 particleDimensions 3 gravity 0 0 -1\n"
        "Species 1\n"
        "LinearViscoelasticSpecies id 0 density 1.9098593 stiffness 2000 dissipation 10\n"
        "Walls 1\n"
        "InfiniteWall id 0 indSpecies 0 position 0 0 0 orientation 0 0 0 1 velocity 0 0 0 angularVelocity 0 0 0 0 force 0 0 0 torque 0 0 0 normal 0 0 -1 factor 1\n"
        "Boundaries 0\n"
        "Particles 0\n"
        "Interactions 0\n"
    );
 
    helpers::writeToFile("Tutorial10.ini.data",
        "1 0 0 0 0 2 2 2\n"
        "1 1 1.5  0 0 0  0.5  0 0 0  0 0 0  0\n"
    );
 
 
    DPMBase Tutorial10;
    //use readRestartFile to load information from a restart file
    Tutorial10.readRestartFile("Tutorial10.ini.restart");
 
    //use readDataFile to load information from a data file
    Tutorial10.readDataFile("Tutorial10.ini.data");
 
    //now start the calculations
    Tutorial10.solve();
 
    return 0;
}

Return to tutorial T10: How to load restart and data files

Axisymmetric walls inside a cylindrical domain (code)

Goto tutorial T11: Axisymmetric walls inside a cylindrical domain

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
//Tutorial 11 introduces the axisymmetrical walls. Two asymmetrical walls are created into a
//cylindrical boundary.
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T9
*/
 
#include "Mercury3D.h"
#include "Walls/IntersectionOfWalls.h"
#include "Walls/AxisymmetricIntersectionOfWalls.h"
#include "Species/LinearViscoelasticSpecies.h"
 
class Tutorial11 : public Mercury3D
{
public:
    Tutorial11(const Mdouble width, const Mdouble height){
        logger(INFO, "Tutorial 11");
        setName("Tutorial11");
        setFileType(FileType::ONE_FILE);
        restartFile.setFileType(FileType::ONE_FILE);
        fStatFile.setFileType(FileType::NO_FILE);
        eneFile.setFileType(FileType::NO_FILE);
        particleHandler.setWriteVTK(true);
        wallHandler.setWriteVTK(FileType::ONE_FILE);
        setXBallsAdditionalArguments("-v0 -solidf");
 
        //specify body forces
        setGravity(Vec3D(0.0, 0.0, -9.8));
 
        //set domain accordingly (domain boundaries are not walls!)
        setXMin(0.0);
        setXMax(width);
        setYMin(0.0);
        setYMax(width);
        setZMin(0.0);
        setZMax(height);
    }
 
    void setupInitialConditions() override {
        Vec3D mid = {
                (getXMin() + getXMax()) / 2.0,
                (getYMin() + getYMax()) / 2.0,
                (getZMin() + getZMax()) / 2.0};
        const Mdouble halfWidth = (getXMax() - getXMin()) / 2.0;
 
        //Cylinder
        AxisymmetricIntersectionOfWalls w1;
        w1.setSpecies(speciesHandler.getObject(0));
        w1.setPosition(Vec3D(mid.X, mid.Y, 0));
        w1.setAxis(Vec3D(0, 0, 1));
        w1.addObject(Vec3D(1, 0, 0), Vec3D((getXMax() - getXMin()) / 2.0, 0, 0));
        // Display the cylinder only, when not coinciding with the neck
        std::vector<Mdouble> displayedSegmentsCylinder = {-constants::inf, mid.Z - contractionHeight, mid.Z + contractionHeight, constants::inf};
        w1.setDisplayedSegments(displayedSegmentsCylinder);
        wallHandler.copyAndAddObject(w1);
 
        //Cone
        AxisymmetricIntersectionOfWalls w2;
        w2.setSpecies(speciesHandler.getObject(0));
        w2.setPosition(Vec3D(mid.X, mid.Y, 0));
        w2.setAxis(Vec3D(0, 0, 1));
        std::vector<Vec3D> points(3);
        //define the neck as a prism through corners of your prismatic wall in clockwise direction
        points[0] = Vec3D(halfWidth, 0.0, mid.Z + contractionHeight);
        points[1] = Vec3D(halfWidth - contractionWidth, 0.0, mid.Z);
        points[2] = Vec3D(halfWidth, 0.0, mid.Z - contractionHeight);
        w2.createOpenPrism(points);
        // Display the neck only until reaching the cylinder
        std::vector<Mdouble> displayedSegmentsPrism = {mid.Z - contractionHeight, mid.Z + contractionHeight};
        w2.setDisplayedSegments(displayedSegmentsPrism);
        wallHandler.copyAndAddObject(w2);
 
        //Flat surface
        InfiniteWall w0;
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(0, 0, -1), Vec3D(0, 0, mid.Z));
        wallHandler.copyAndAddObject(w0);
 
        const int N = 600; //maximum number of particles
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        for (Mdouble z = mid.Z + contractionHeight;
             particleHandler.getNumberOfObjects() <= N;
             z += 2.0 * maxParticleRadius)
        {
            for (Mdouble r = halfWidth - maxParticleRadius; r > 0; r -= 1.999 * maxParticleRadius)
            {
                for (Mdouble c = 2.0 * maxParticleRadius; c <= 2 * constants::pi * r; c += 2.0 * maxParticleRadius)
                {
                    p0.setRadius(random.getRandomNumber(minParticleRadius, maxParticleRadius));
                    p0.setPosition(Vec3D(mid.X + r * mathsFunc::sin(c / r),
                                         mid.Y + r * mathsFunc::cos(c / r),
                                         z + p0.getRadius()));
                    const Mdouble vz = random.getRandomNumber(-1, 0);
                    const Mdouble vx = random.getRandomNumber(-1, 1);
                    p0.setVelocity(Vec3D(vx,0.0,vz));
                    particleHandler.copyAndAddObject(p0);
                }
            }
        }
    }
 
    //Initially, a wall is inserted in the neck of the hourglass to prevent particles flowing through.
    //This wall is moved to form the base of the hourglass at time 0.9
    void actionsAfterTimeStep() override {
        if (getTime() < 0.9 && getTime() + getTimeStep() > 0.9)
        {
            logger(INFO, "Shifting bottom wall downward");
            dynamic_cast<InfiniteWall*>(wallHandler.getLastObject())->set(Vec3D(0, 0, -1), Vec3D(0, 0, getZMin()));
        }
    }
 
    Mdouble contractionWidth{};
    Mdouble contractionHeight{};
    Mdouble minParticleRadius{};
    Mdouble maxParticleRadius{};
 
};
 
int main(int argc, char *argv[])
{
    Mdouble width = 10e-2; // 10cm
    Mdouble height = 60e-2; // 60cm
 
    //Point the object HG to class
    Tutorial11 HG(width,height);
 
    //Specify particle radius:
    //these parameters are needed in setupInitialConditions()
    HG.minParticleRadius = 6e-3; // 6mm
    HG.maxParticleRadius = 10e-3; //10mm
 
    //Specify the number of particles
 
    //specify how big the wedge of the contraction should be
    const Mdouble contractionWidth = 2.5e-2; //2.5cm
    const Mdouble contractionHeight = 5e-2; //5cm
    HG.contractionWidth = contractionWidth;
    HG.contractionHeight = contractionHeight;
 
    //make the species of the particle and wall
    LinearViscoelasticSpecies species;
    species.setDensity(2000);
 
    species.setStiffness(1e5);
    species.setDissipation(0.63);
    HG.speciesHandler.copyAndAddObject(species);
 
 
    //test normal forces
    const Mdouble minParticleMass = species.getDensity() * 4.0 / 3.0 * constants::pi * mathsFunc::cubic(HG.minParticleRadius);
    //Calculates collision time for two copies of a particle of given dissipation_, k, effective mass
    logger(INFO, "minParticleMass = %", minParticleMass);
    //Calculates collision time for two copies of a particle of given dissipation_, k, effective mass
    const Mdouble tc = species.getCollisionTime(minParticleMass);
    logger(INFO, "tc  = %", tc);
    //Calculates restitution coefficient for two copies of given dissipation_, k, effective mass
    const Mdouble r = species.getRestitutionCoefficient(minParticleMass);
    logger(INFO, "restitution coefficient = %", r);
 
    //time integration parameters
    HG.setTimeStep(tc / 10.0);
    HG.setTimeMax(5.0);
    HG.setSaveCount(500);
 
    HG.solve(argc, argv);
    return 0;
}

Return to tutorial T11: Axisymmetric walls inside a cylindrical domain

Prismatic walls using points (code)

Goto tutorial T12: Creating objects by using Prisms

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
// Tutorial 12
 
/*
** This file is annotated with DoxyFile comments in order to show the code on
** the documentation - This is not needed for your real drivers.
** Please ignore these comments.
**
** For full documentation of this code, go to http://docs.mercurydpm.org/Alpha/d0/db0/BeginnerTutorials.html#T8
*/
 
#include <Mercury3D.h>
#include <Walls/InfiniteWall.h>
#include <Walls/IntersectionOfWalls.h>
#include <Species/LinearViscoelasticSpecies.h>
 
class Tutorial12 : public Mercury3D
{
public:
    void setupInitialConditions() override {
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setRadius(0.005);
        p0.setPosition(Vec3D(0.15 * getXMax(), 0.3335 * getYMax(), 0.0));
        p0.setVelocity(Vec3D(1.2, 0.2, 0.0));
        particleHandler.copyAndAddObject(p0);
        
        // Creating outer walls
        InfiniteWall w0;
        
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(1.0, 0.0, 0.0), Vec3D(getXMax(), 0.0, 0.0));
        wallHandler.copyAndAddObject(w0);
        
        w0.set(Vec3D(-1.0, 0.0, 0.0), Vec3D(getXMin(), 0.0, 0.0));
        wallHandler.copyAndAddObject(w0);
        
        w0.set(Vec3D(0.0, 1.0, 0.0), Vec3D(0.0, getYMax(), 0.0));
        wallHandler.copyAndAddObject(w0);
        
        w0.set(Vec3D(0.0, -1.0, 0.0), Vec3D(0.0, getYMin(), 0.0));
        wallHandler.copyAndAddObject(w0);
 
        // Defining the object in het center of the box
        // Create an intersection of walls object w1
        IntersectionOfWalls w1;
        // Set the species of the wall
        w1.setSpecies(speciesHandler.getObject(0));
        std::vector<Vec3D> Points(4);
        // Define the vertices of the insert and ensure that points are in clockwise order
        Points[0] = Vec3D(0.25 * getXMax(), -0.25 * getYMax(), 0.0);
        Points[1] = Vec3D(-0.25 * getXMax(), -0.25 * getYMax(), 0.0);
        Points[2] = Vec3D(-0.25 * getXMax(), 0.25 * getYMax(), 0.0);
        Points[3] = Vec3D(0.25 * getXMax(), 0.25 * getYMax(), 0.0);
        // Creating the object from the vector containing points
        w1.createPrism(Points);
        /* Define the angular velocity of the object (optional).
         * Setting the angular velocity of the object results in rotation of the object around
         * the origin (0.0,0.0,0.0). If the object is build such that the center of rotation of the object
         * coincides with the origin the object will rotate around its own axis. */
        w1.setAngularVelocity(Vec3D(0.0,0.0,1.0));
        /* Define the translational velocity of the object (optional) */
        w1.setVelocity(Vec3D(0.0,0.0,0.0));
        /* Set the desired position of center of the object (optional).
         * With this command you can position your object (with the set angular velocity)
         * at a specific location in the domain.*/
        w1.setPosition(Vec3D(0.5 * getXMax(),0.5 * getYMax(),0.0));
        // Add the object to wall w1
        wallHandler.copyAndAddObject(w1);
        
    }
};
 
int main(int argc, char* argv[])
{
    // Problem setup
    Tutorial12 problem; // instantiate an object of class Tutorial 12
    
    problem.setName("Tutorial12");
    problem.setGravity(Vec3D(0.0, 0.0, 0.0));
    problem.setXMax(0.5);
    problem.setYMax(0.5);
    problem.setZMax(0.1);
    problem.setTimeMax(5.0);
 
    // The normal spring stiffness and normal dissipation is computed and set as
    // For collision time tc=0.005 and restitution coefficient rc=1.0,
    LinearViscoelasticSpecies species;
    species.setDensity(2500.0); //sets the species type_0 density
    species.setStiffness(258.5);//sets the spring stiffness.
    species.setDissipation(0.0); //sets the dissipation.
    problem.speciesHandler.copyAndAddObject(species);
 
    
    problem.setSaveCount(10);
    problem.dataFile.setFileType(FileType::ONE_FILE);
    problem.restartFile.setFileType(FileType::ONE_FILE);
    problem.fStatFile.setFileType(FileType::NO_FILE);
    problem.eneFile.setFileType(FileType::NO_FILE);
    problem.particleHandler.setWriteVTK(true);
    problem.wallHandler.setWriteVTK(FileType::ONE_FILE);
 
    problem.setTimeStep(.005 / 50.0); // (collision time)/50.0
    
    problem.solve(argc, argv);
    
    return 0;
}

Return to tutorial T12: Creating objects by using Prisms

Particles driven by a rotating coil (code)

Goto tutorial Particles driven by a rotating coil

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#include "Mercury3D.h"
#include "Walls/InfiniteWall.h"
#include "Walls/InfiniteWallWithHole.h"
#include "Walls/Coil.h"
#include "Species/LinearViscoelasticSpecies.h"
 
 
class CoilSelfTest : public Mercury3D
{
private:
    void setupInitialConditions() override
    {
        InfiniteWall w;
        w.setSpecies(speciesHandler.getObject(0));
        //front wall
        w.set(Vec3D(-1, 0, 0), Vec3D(getXMin(), 0, 0));
        wallHandler.copyAndAddObject(w);
 
        //back wall
        w.set(Vec3D(1, 0, 0), Vec3D(getXMax(), 0, 0));
        wallHandler.copyAndAddObject(w);
 
        //bottom wall
        w.set(Vec3D(0, -1, 0), Vec3D(0, getYMin(), 0));
        wallHandler.copyAndAddObject(w);
 
        //top wall
        w.set(Vec3D(0, 1, 0), Vec3D(0, getYMax(), 0));
        wallHandler.copyAndAddObject(w);
 
        //left wall
        w.set(Vec3D(0, 0, -1), Vec3D(0, 0, getZMin()));
        wallHandler.copyAndAddObject(w);
 
        //right wall
        InfiniteWallWithHole rightWall;
        rightWall.setSpecies(speciesHandler.getObject(0));
        rightWall.set(Vec3D(0, 0, 1), getZMax(), 1.0);
        wallHandler.copyAndAddObject(rightWall);
 
        // creation of the coil and setting of its properties
        Coil coil;
        coil.setSpecies(speciesHandler.getObject(0));
        // the syntax to set the coil geometry is as follows:
        // set(Start position, Length, Radius, Number of turns, Rotation speed, Thickness)
        coil.set(Vec3D(0, 0, 0), 1.0, 1.0 - particleRadius, 2.0, -1.0, 0.5 * particleRadius);
        auto pCoil = wallHandler.copyAndAddObject(coil);
 
        particleHandler.clear();
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
        p0.setRadius(particleRadius);
        /*
        //Single test case
        double distance;
                Vec3D normal;
        p0.setPosition(Vec3D(1.0,0.0,0.0));
        if(coil->getDistance_and_normal(p0, distance, normal))
                        std::cout<<"Collision, distance screw="<<distance<<std::endl;
                else
                        std::cout<<"No collision, distance screw="<<distance<<std::endl;
         */
 
        // Nx*Ny*Nz particles are created and placed on evenly spaced positions in
        // the domain [xMin_,xMax_]*[yMin_,yMax_]*[zMin_,zMax_] (unless their position
        // is already occupied by the coil).
        const unsigned int Nx = static_cast<unsigned int> (std::floor((getXMax() - getXMin()) / (2.1 * particleRadius)));
        const unsigned int Ny = static_cast<unsigned int> (std::floor((getYMax() - getYMin()) / (2.1 * particleRadius)));
        const unsigned int Nz = static_cast<unsigned int> (std::floor((getZMax() - getZMin()) / (2.1 * particleRadius)));
        Mdouble distance;
        Vec3D normal;
 
        for (unsigned int i = 0; i < Nx; i++)
        {
            for (unsigned int j = 0; j < Ny; j++)
            {
                for (unsigned int k = 0; k < Nz; k++)
                {
                    p0.setPosition(Vec3D(getXMin() + (getXMax() - getXMin()) * (0.5 + i) / Nx,
                                         getYMin() + (getYMax() - getYMin()) * (0.5 + j) / Ny,
                                         getZMin() + (getZMax() - getZMin()) * (0.5 + k) / Nz));
                    if (!pCoil->getDistanceAndNormal(p0, distance, normal)) //if there is no collision with the coil
                    {
                        particleHandler.copyAndAddObject(p0);
                    }
                    else
                    {
                        logger(DEBUG, "particle at position % could not be inserted", p0.getPosition());
                    }
                }
            }
        }
    }
    void actionsBeforeTimeStep() override
    {
        if (getTime() > 1)
        {
            coil->move_time(getTimeStep());
        }
    }
 
public:
    // ! [CST:datamembers]
    Coil* coil;
    Mdouble particleRadius;
};
 
int main(int argc UNUSED, char* argv[] UNUSED)
{
    
    // create CoilSelfTest object
    CoilSelfTest problem;
 
    // set some basic problem properties
    problem.setName("CoilSelfTest");
    problem.setSystemDimensions(3);
    problem.setGravity(Vec3D(0.0, -9.8, 0.0));
 
    // set problem geometry
    problem.setXMax(1.0);
    problem.setYMax(5.0);
    problem.setZMax(2.0);
    problem.setXMin(-1.0);
    problem.setYMin(-1.0);
    problem.setTimeMax(0.5);
    problem.particleRadius = 0.2;
 
    LinearViscoelasticSpecies species;
    species.setDensity(1000);
    Mdouble tc = 0.05;
    Mdouble restitutionCoefficient = 0.8;
 
    Mdouble particleMass = pow(problem.particleRadius, 3) * constants::pi * 4.0 / 3.0 * species.getDensity();
    species.setCollisionTimeAndRestitutionCoefficient(tc, restitutionCoefficient, particleMass);
    problem.speciesHandler.copyAndAddObject(species);
 
    problem.setTimeStep(0.02 * 0.05);
    problem.setSaveCount(helpers::getSaveCountFromNumberOfSavesAndTimeMaxAndTimeStep(1000, problem.getTimeMax(),
                                                                                     problem.getTimeStep()));
    problem.solve();
    
}

Return to tutorial Particles driven by a rotating coil

Particles on an inclined chute (code)

Goto tutorial Particles on an inclined chute

#include <iostream>
#include <Species/LinearViscoelasticSpecies.h>
#include "Chute.h"
 
// Creates a quasi-2D inclined plane with inflow conditions on the left boundary, 
// and deletion of particles when they exit the domain on the right.
int main() 
{
    
    // Problem parameters
    Chute problem;
    
    problem.setName("ChuteDemo");   // data output file name
    problem.setGravity({0,-1,0});
    problem.setSaveCount(102);      // number of time steps skipped between saves
    Mdouble tc = 2.5e-3;            // collision time
    problem.setTimeStep(0.02 * tc); // actual time step
    problem.setTimeMax(0.5);        // maximum time
                                    // NB: number of time steps saved to output files
                                    // is timeMax/(timeStep*saveCount)
    
    // Particle radii
    problem.setFixedParticleRadius(0.001);  // radius of fixed chute bottom particles
    problem.setInflowParticleRadius(0.001); // radius of (monodisperse) inflow particles
 
    // Particle species
    LinearViscoelasticSpecies species;              // initialise species
    species.setHandler(&problem.speciesHandler);    // assign problem species handler to species
    species.setDensity(2000);                       // particle density
    species.setCollisionTimeAndRestitutionCoefficient( tc, 
            0.8, species.getMassFromRadius(problem.getInflowParticleRadius())); // material properties
    problem.speciesHandler.copyAndAddObject(species);   // assign species to problem species handler
    
 
    // Chute properties
    problem.setChuteAngle(30.0);                    // set angle of chute relative to horizontal
    problem.setXMax(0.1);                           // chute length = 0.1
    problem.setYMax(2.0 * problem.getInflowParticleRadius()); // chute width = 1 particle diameter
 
    // Inflow properties
    problem.setInflowHeight(0.1);                   // particle inflow between 0 <= Z <= 0.1
    problem.setInflowVelocity(0.1);                 // particle inflow mean velocity
    problem.setInflowVelocityVariance(0.02);        // particle inflow velocity variance (in ratio of the mean velocity)
    //Write paraview data
    problem.particleHandler.setWriteVTK(true);
    problem.wallHandler.setWriteVTK(true);
 
 
    /*problem.particleHandler.setWriteVTK(true);
    problem.wallHandler.setWriteVTK(true);*/
 
 
    //solve
    problem.solve();
} // the end

Return to tutorial Particles on an inclined chute

Particles on a chute with a multilayered bottom (code)

Goto tutorial Particles on a chute with a multilayered bottom

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#include <sstream>
#include <Species/LinearViscoelasticSpecies.h>
#include "Chute.h"
 
int main()
{
    //Print description
    logger(INFO, "\nDemo of the chute flow with a rough bottom");
    
    //Construct the problem and assign a name
    Chute roughBottomSelfTest;
    roughBottomSelfTest.setName("roughBottomMultiLayer");
 
 
    //Set time stepping parameters
    roughBottomSelfTest.setTimeStep(1e-4);
    roughBottomSelfTest.setTimeMax(0.1);
 
    //Set output parameter: write to the output files every 50 time steps
    roughBottomSelfTest.setSaveCount(50);
 
    //Set the particle radii and the type of the rough bottom
    roughBottomSelfTest.setInflowParticleRadius(0.5);
    roughBottomSelfTest.setFixedParticleRadius(0.5);
    roughBottomSelfTest.setRoughBottomType(MULTILAYER);
 
    //Set the particle properties
    LinearViscoelasticSpecies* species = roughBottomSelfTest.speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
    species->setDensity(6.0/constants::pi);
    species->setStiffness(2e5);
    species->setDissipation(25.0);
 
    //Set the chute properties
    roughBottomSelfTest.setChuteAngleAndMagnitudeOfGravity(24.0, 1.0);
    roughBottomSelfTest.setChuteLength(3.0);
    roughBottomSelfTest.setChuteWidth(3.0);
 
    //Solve the system
    roughBottomSelfTest.solve();
}

Return to tutorial Particles on a chute with a multilayered bottom

Isotropic compression of a cuboidal REV (code)

Goto tutorial Isotropic compression of a cuboidal REV

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
//#include <Species/Species.h>
#include <Mercury3D.h>
#include <Species/LinearViscoelasticSpecies.h>
#include <Boundaries/StressStrainControlBoundary.h>
#include <Boundaries/CubeInsertionBoundary.h>
 
class StressStrainControl : public Mercury3D
{
public:
    //Create the class and passing the user inputs using the constructors
    StressStrainControl(const Matrix3D& stressGoal, const Matrix3D& strainRate, const Matrix3D& gainFactor,
                        bool isStrainRateControlled)
            : stressGoal_(stressGoal), strainRate_(strainRate), gainFactor_(gainFactor),
              isStrainRateControlled_(isStrainRateControlled)
    {
        //Define the domain size
        setMin(Vec3D(0, 0, 0));
        setMax(Vec3D(1, 1, 1));
        //Calculate the mass from smallest particle
        double mass = 2000 * constants::pi * mathsFunc::cubic(2.0 * 0.08) / 6.0;
        //Define the species for the particles
        LinearViscoelasticSpecies species;
        species.setDensity(2000);
        species.setStiffnessAndRestitutionCoefficient(10000, 0.4, mass);
        speciesHandler.copyAndAddObject(species);
    }
 
private:
    //Here we set up the system parameters, adding particles and define boundaries
    void setupInitialConditions() override
    {
        //Set up the micro parameters.
        double rhop = 2000;     //particle density
        double K1 = 10000;      //particle stiffness
        double en = 0.4;        //restitution coefficient
        double radius = 0.08;   //particle radius
        //Calculate the mass from smallest particle
        double mass = rhop * constants::pi * mathsFunc::cubic(2.0 * radius) / 6.0;
        //Calculate the contact duration
        double tc = std::sqrt(mass / 2.0 / K1 * (mathsFunc::square(constants::pi) +
                                                 mathsFunc::square(log(en))));
        setSystemDimensions(3); //set the 3d simulation
        setTimeStep(tc / 50); //set the timestep based on the shortest contact duration
        setGravity(Vec3D(0.0, 0.0, 0.0)); //set gravity in the system
 
        //Add particles
        CubeInsertionBoundary insertion;
        SphericalParticle particle(speciesHandler.getObject(0));
        particle.setRadius(radius);
        //Insert 96 particles in the subvolume = 0.8*0.8*0.8 = 0.512, if failed by 100 times, it stops
        insertion.set(&particle, 100, 0.8 * getMin(), 0.8 * getMax(), Vec3D(0, 0, 0), Vec3D(0, 0, 0));
        insertion.setInitialVolume(1.0);
        insertion.checkBoundaryBeforeTimeStep(this);
        logger(INFO,"Inserted % particles",particleHandler.getSize());
 
        //Delete all existing boundaries
        boundaryHandler.clear();
        //Set up the deformation mode using the boundaries that are based on user's input
        StressStrainControlBoundary boundary;
        boundary.setHandler(&boundaryHandler);
        boundary.set(stressGoal_, strainRate_, gainFactor_, isStrainRateControlled_);
        boundaryHandler.copyAndAddObject(boundary);
 
        
    }
 
    //initialize the private variables for passing the user's inputs in the constructor
    Matrix3D stressGoal_;
    Matrix3D strainRate_;
    Matrix3D gainFactor_;
    bool isStrainRateControlled_;
};
 
//Here we define all the control parameters and solve the problem
int main(int argc UNUSED, char* argv[] UNUSED)
{
    //We want strainrate control, so here we set the target Stress to free, all to zero
    Matrix3D stressGoal;
    stressGoal.XX = 0.0;
    stressGoal.YY = 0.0;
    stressGoal.ZZ = 0.0;
    stressGoal.XY = 0.0;
 
    //Here we set the isotropic strainrate tensor, negative sign means compression
    Matrix3D strainRate;
    strainRate.XX = -0.02;
    strainRate.YY = -0.02;
    strainRate.ZZ = -0.02;
    strainRate.XY = 0.0;
 
    //This is default gainFactor, if you choose to use stress control, you might want to tune this one
    Matrix3D gainFactor;
    gainFactor.XY = 0.0001;
    gainFactor.XX = 0.0001;
    gainFactor.YY = 0.0001;
    gainFactor.ZZ = 0.0001;
 
    //This is by default set to true, so particles are controlled by affine movement.
    //If set it to false, then the particles will only follow the boundary movement.
    bool isStrainRateControlled = true;
    
    //Define the object and problem to solve
    StressStrainControl dpm(stressGoal, strainRate, gainFactor, isStrainRateControlled);
    
    //Set name
    dpm.setName("Tutorial_IsotropicCompression");
    //Set output and time stepping properties
    dpm.setTimeMax(10);
    //Set the SaveCount, i.e. 100 timesteps saving one snapshot
    dpm.setSaveCount(100);
    //Currently all the normal file outputs are switched off, simply switch it on by replacing NO_FILE to ONE_FILE
//    dpm.dataFile.setFileType(FileType::NO_FILE);
//    dpm.restartFile.setFileType(FileType::NO_FILE);
//    dpm.fStatFile.setFileType(FileType::NO_FILE);
//    dpm.eneFile.setFileType(FileType::NO_FILE);
    //Set output the particle information in VTK for ParaView Visualisation
    dpm.particleHandler.setWriteVTK(true);
    //Because of periodic boundary, out put wall files is not necessary in this case
    //dpm.wallHandler.setWriteVTK(true);
    //Solve the problem
    dpm.solve();
}

Return to tutorial Isotropic compression of a cuboidal REV

Pure shear (constant volume) of a cuboidal REV (code)

Goto tutorial Pure shear (constant volume) of a cuboidal REV

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
//#include <Species/Species.h>
#include <Mercury3D.h>
#include <Species/LinearViscoelasticSpecies.h>
#include <Boundaries/StressStrainControlBoundary.h>
#include <Boundaries/CubeInsertionBoundary.h>
 
class StressStrainControl : public Mercury3D
{
public:
    //Create the class and passing the user inputs using the constructors
    StressStrainControl(const Matrix3D& stressGoal, const Matrix3D& strainRate, const Matrix3D& gainFactor,
                        bool isStrainRateControlled)
            : stressGoal_(stressGoal), strainRate_(strainRate), gainFactor_(gainFactor),
              isStrainRateControlled_(isStrainRateControlled)
    {
        //Define the domain size
        setMin(Vec3D(0, 0, 0));
        setMax(Vec3D(1, 1, 1));
        //Calculate the mass from smallest particle
        double mass = 2000 * constants::pi * mathsFunc::cubic(2.0 * 0.08) / 6.0;
        //Define the species for the particles
        LinearViscoelasticSpecies species;
        species.setDensity(2000);
        species.setStiffnessAndRestitutionCoefficient(10000, 0.4, mass);
        speciesHandler.copyAndAddObject(species);
    }
 
private:
    //Here we set up the system parameters, adding particles and define boundaries
    void setupInitialConditions() override
    {
        //Set up the micro parameters.
        double rhop = 2000;     //particle density
        double K1 = 10000;      //particle stiffness
        double en = 0.4;        //restitution coefficient
        double radius = 0.05;   //particle radius
        //Calculate the mass from smallest particle
        double mass = rhop * constants::pi * mathsFunc::cubic(2.0 * radius) / 6.0;
        //Calculate the contact duration
        double tc = std::sqrt(mass / 2.0 / K1 * (mathsFunc::square(constants::pi) +
                                                 mathsFunc::square(log(en))));
        setSystemDimensions(3); //set the 3d simulation
        setTimeStep(tc / 50); //set the timestep based on the shortest contact duration
        setGravity(Vec3D(0.0, 0.0, 0.0)); //set gravity in the system
 
        //Add particles
        CubeInsertionBoundary insertion;
        SphericalParticle particle(speciesHandler.getObject(0));
        //Insert 96 particles in the subvolume = 1.0*1.0*1.0 = 1.0, if failed by 100 times, it stops
        insertion.set(&particle, 100, getMin(), getMax(), Vec3D(0, 0, 0), Vec3D(0, 0, 0));
        insertion.setInitialVolume(1.0);
        insertion.checkBoundaryBeforeTimeStep(this);
        logger(INFO,"Inserted % particles",particleHandler.getSize());
 
        //Set up the deformation mode using the boundaries that are based on user's input
        boundaryHandler.clear();                  // Delete all exist boundaries
        StressStrainControlBoundary boundary;
        boundary.setHandler(&boundaryHandler);
        boundary.set(stressGoal_, strainRate_, gainFactor_, isStrainRateControlled_);
        boundaryHandler.copyAndAddObject(boundary);
 
        
    }
 
    //Initialize the private variables for passing the user's inputs in the constructor
    Matrix3D stressGoal_;
    Matrix3D strainRate_;
    Matrix3D gainFactor_;
    bool isStrainRateControlled_;
};
 
//Here we define all the control parameters and solve the problem
int main(int argc UNUSED, char* argv[] UNUSED)
{
    //We want strainrate control, so here we set the target Stress to free, all to zero
    Matrix3D stressGoal;
    stressGoal.XX = 0.0;
    stressGoal.YY = 0.0;
    stressGoal.ZZ = 0.0;
    stressGoal.XY = 0.0;
 
    //Here we set the pure shear strainrate tensor, negative sign means compression
    Matrix3D strainRate;
    strainRate.XX = 0.4;
    strainRate.YY = -0.4;
    strainRate.ZZ = 0.0;
    strainRate.XY = 0.0;
 
    //This is default gainFactor, if you choose to use stress control, you might want to tune this one
    Matrix3D gainFactor;
    gainFactor.XY = 0.0001;
    gainFactor.XX = 0.0001;
    gainFactor.YY = 0.0001;
    gainFactor.ZZ = 0.0001;
 
    //This is by default set to true, so particles are controlled by affine movement.
    //If set it to false, then the particles will only follow the boundary movement.
    bool isStrainRateControlled = true;
    
    //Define the object and problem to solve
    StressStrainControl dpm(stressGoal, strainRate, gainFactor, isStrainRateControlled);
    
    //Set name
    dpm.setName("Tutorial_PureShear");
    //Set output and time stepping properties
    dpm.setTimeMax(1);
    //Set the SaveCount, i.e. 100 timesteps saving one snapshot
    dpm.setSaveCount(10);
    //Currently all the normal file outputs are switched off, simply switch it on by replacing NO_FILE to ONE_FILE
    dpm.dataFile.setFileType(FileType::NO_FILE);
    dpm.restartFile.setFileType(FileType::NO_FILE);
    dpm.fStatFile.setFileType(FileType::NO_FILE);
    dpm.eneFile.setFileType(FileType::NO_FILE);
    //Set output the particle information in VTK for ParaView Visualisation
    dpm.particleHandler.setWriteVTK(true);
    //Because of periodic boundary, out put wall files is not necessary in this case
    //dpm.wallHandler.setWriteVTK(true);
    //Solve the problem
    dpm.solve();
}

Return to tutorial Pure shear (constant volume) of a cuboidal REV

(Rigid_clumps) (code)

Goto tutorial (Rigid_clumps) (code)

Return to tutorial (Rigid_clumps) (code)

or Dzhanibekov theorem (code)

Goto tutorial Tennisracket or Dzhanibekov theorem

Return to tutorial Tennisracket or Dzhanibekov theorem

effect (code)

Goto tutorial Domino effect

Return to tutorial Domino effect

(code)

Goto tutorial Simple shear (constant stress) of a cuboidal REV

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
//#include <Species/Species.h>
#include <Mercury3D.h>
#include <Species/LinearViscoelasticSpecies.h>
#include <Boundaries/StressStrainControlBoundary.h>
#include <Boundaries/CubeInsertionBoundary.h>
 
class StressStrainControl : public Mercury3D
{
public:
    //Create the class and passing the user inputs using the constructors
    StressStrainControl(const Matrix3D& stressGoal, const Matrix3D& strainRate, const Matrix3D& gainFactor,
                        bool isStrainRateControlled)
            : stressGoal_(stressGoal), strainRate_(strainRate), gainFactor_(gainFactor),
              isStrainRateControlled_(isStrainRateControlled)
    {
        //Define the domain size
        setMin(Vec3D(0, 0, 0));
        setMax(Vec3D(1, 1, 1));
        //Calculate the mass from smallest particle
        double mass = 2000 * constants::pi * mathsFunc::cubic(2.0 * 0.05) / 6.0;
        //Define the species for the particles
        LinearViscoelasticSpecies species;
        species.setDensity(2000);
        species.setStiffnessAndRestitutionCoefficient(10000, 0.4, mass);
        speciesHandler.copyAndAddObject(species);
    }
 
private:
    //Here we set up the system parameters, adding particles and define boundaries
    void setupInitialConditions() override
    {
        //Set up the micro parameters.
        double rhop = 2000;     //particle density
        double K1 = 10000;      //particle stiffness
        double en = 0.4;        //restitution coefficient
        double radius = 0.05;   //particle radius
        //Calculate the mass from smallest particle
        double mass = rhop * constants::pi * mathsFunc::cubic(2.0 * radius) / 6.0;
        //Calculate the contact duration
        double tc = std::sqrt(mass / 2.0 / K1 * (mathsFunc::square(constants::pi) +
                                                 mathsFunc::square(log(en))));
        setSystemDimensions(3); //set the 3d simulation
        setTimeStep(tc / 50); //set the timestep based on the shortest contact duration
        setGravity(Vec3D(0.0, 0.0, 0.0)); //set gravity in the system
 
        //Add particles
        CubeInsertionBoundary insertion;
        SphericalParticle particle(speciesHandler.getObject(0));
        //Insert 96 particles in the subvolume = 1.0*1.0*1.0 = 1.0, if failed by 100 times, it stops
        insertion.set(&particle, 100, getMin(), getMax(), Vec3D(0, 0, 0), Vec3D(0, 0, 0));
        insertion.setInitialVolume(1.0);
        insertion.checkBoundaryBeforeTimeStep(this);
        logger(INFO,"Inserted % particles",particleHandler.getSize());
 
        //Set up the deformation mode using the boundaries that are based on user's input
        boundaryHandler.clear();                  // Delete all exist boundaries
        StressStrainControlBoundary boundary;
        boundary.setHandler(&boundaryHandler);
        boundary.set(stressGoal_, strainRate_, gainFactor_, isStrainRateControlled_);
        boundaryHandler.copyAndAddObject(boundary);
 
        
    }
 
    //Initialize the private variables for passing the user's inputs in the constructor
    Matrix3D stressGoal_;
    Matrix3D strainRate_;
    Matrix3D gainFactor_;
    bool isStrainRateControlled_;
};
 
//Here we define all the control parameters and solve the problem
int main(int argc UNUSED, char* argv[] UNUSED)
{
    //Here we want to set stress in y direction to constant, 100 is just an example,
    //if you would like to keep the volume constant, then everything here should be set to zero
    Matrix3D stressGoal;
    stressGoal.XX = 0.0;
    stressGoal.YY = 100.0;
    stressGoal.ZZ = 0.0;
    stressGoal.XY = 0.0;
 
    //Here we set the simple shear strainrate tensor, negative sign means compression
    Matrix3D strainRate;
    strainRate.XX = 0.0;
    strainRate.YY = 0.0;
    strainRate.ZZ = 0.0;
    strainRate.XY = 1.0;
 
    //This is default gainFactor, if you choose to use stress control, you might want to tune this one
    Matrix3D gainFactor;
    gainFactor.XY = 0.0001;
    gainFactor.XX = 0.0001;
    gainFactor.YY = 0.0001;
    gainFactor.ZZ = 0.0001;
 
    //This is by default set to true, so particles are controlled by affine movement.
    //If set it to false, then the particles will only follow the boundary movement.
    bool isStrainRateControlled = true;
    
    //Define the object and problem to solve
    StressStrainControl dpm(stressGoal, strainRate, gainFactor, isStrainRateControlled);
    
    //Set name
    dpm.setName("Tutorial_SimpleShear");
    //Set output and time stepping properties
    dpm.setTimeMax(10);
    //Set the SaveCount, i.e. 100 timesteps saving one snapshot
    dpm.setSaveCount(100);
    //Currently all the normal file outputs are switched off, simply switch it on by replacing NO_FILE to ONE_FILE
    dpm.dataFile.setFileType(FileType::NO_FILE);
    dpm.restartFile.setFileType(FileType::NO_FILE);
    dpm.fStatFile.setFileType(FileType::NO_FILE);
    dpm.eneFile.setFileType(FileType::NO_FILE);
    //Set output the particle information in VTK for ParaView Visualisation
    dpm.particleHandler.setWriteVTK(true);
    //Because of periodic boundary, out put wall files is not necessary in this case
    //dpm.wallHandler.setWriteVTK(true);
    //Solve the problem
    dpm.solve();
}

Return to tutorial Simple shear (constant stress) of a cuboidal REV