Beginner tutorials

This page contains the following tutorials:

• T1: Particle motion in outer space
• T2: Particle motion on earth
• T3: Bouncing ball (elastic)
• T4: Bouncing ball with dissipation (inelastic)
• T5: Elastic collision (2 particles)
• T6: Elastic collisions with periodic boundaries
• T7: Motion of a particle in a two dimensional (2D) box
• T8: Motion of a particle in a box with an obstacle
• T9: Motion of a ball over an inclined plane (Sliding + Rolling)
• T11: Axisymmetric walls inside a cylindrical domain
• T12: Creating objects by using Prisms

T1: Particle motion in outer space

Particle moving with a constant velocity in outer space.

Problem description:

The first tutorial is setup to simulate a particle moving with a constant velocity in the absence of gravity. You can find the simulation setup in Tutorial1_ParticleInOuterSpace.cpp The detailed description of this tutorial is presented below.

Headers:

The headers are necessary to set up the simulations. They come from the kernel and standard libraries (see Developing a kernel feature). Headers files usually have a .h or .hpp extensions. For example:

#include <Mercury3D.h>
#include <Species/Species.h>
#include <Species/LinearViscoelasticSpecies.h>

Before main function:

A class is the guide for objects to connect with handlers.
You can find in the code:

class Tutorial1 : public Mercury3D 

Here is defined the class Tutorial1, which is inherited from the main class Mercury3D. Then, you can call the function to fit the initial conditions of the problem:

void setupInitialConditions()

It creates the problems’ initial conditions to define the radius, position, velocity of particles. In order to create a particle, the following structure must be used:

        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setRadius(0.05); // sets particle radius
        p0.setPosition(Vec3D(0.1 * getXMax(), 0.1 * getYMax(), 0.1 * getZMax())); // sets particle position
        p0.setVelocity(Vec3D(0.5, 0.1, 0.1));// sets particle velocity
        particleHandler.copyAndAddObject(p0);

The particle data members as positions, radius and velocities are pointed by using an object (P0) to their class member. For the case of the mentioned data members, the class member is called BaseParticle. Then, the defined information is sent to its corresponding handler: particleHandler by using the function copyAndAddObject

Main function:

In the main function, the global parameters of the problem are defined. It includes gravity, spatial dimensions (x,y,z), total run time, the type of contact law, etc.
The object 'problem' is an instance of the defined class Tutorial1. Global parameters and species will point to the object 'problem'.

Tutorial1 problem;

Then, the data members are defined as:

    problem.setName("Tutorial1");
    problem.setSystemDimensions(3);
    problem.setGravity(Vec3D(0.0, 0.0, 0.0));
    problem.setXMax(1.0);
    problem.setYMax(1.0);
    problem.setZMax(1.0);
    problem.setTimeMax(2.0);

Subsequently, the Species of the problem are defined. It means, the properties (density,stiffness) and the corresponding contact law (for this tutorial, see LinearViscoelasticSpecies). Initially, when a particle is created, it attains the properties of a default species type with ‘0’ as its index.

    //Set the species of particles and walls
    //The normal spring stiffness and normal dissipation is computed and set as
    //For collision time tc=0.005 and restitution coefficeint 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);

Here, 'species' is an object of the class: LinearViscoelasticSpecies. Data members are fitted to point to this object, and then, coppied and added to their corresponding handler speciesHandler.

Outputs:
Data output is vital to analyse simulations, which leads to defining ways to save the simulation data for post processing.
The simulations generate several types of data files.

    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);
    logger(INFO, "run number: %", problem.dataFile.getCounter());

Thereby, data members and member functions are pointed to the object 'problem'.
For XBalls users, additional display options can be set as below

    problem.setXBallsAdditionalArguments("-solidf -v0");

After all the simulation parameters are set, we reach a point where we put all the above bits of code in action by the following statements

    problem.setTimeStep(0.005 / 50.0); // (collision time)/50.0
    problem.solve(argc, argv);

T2: Particle motion on earth

Particle falling due to gravity.

Problem description:

In Tutorial2_ParticleUnderGravity.cpp, we simulate a particle when dropped under the influence of gravity, \( g = 9.81 m/s^2\). Basically, this tutorial is an extension of Tutorial1_ParticleInOuterSpace.cpp with few minor changes.
All we need to do is to change the initial particle position and velocity in the inherited class Tutorial2.

        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);

Also, the gravity vector in the main function:

    problem.setName("Tutorial2");
    problem.setSystemDimensions(3);
    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);

T3: Bouncing ball (elastic)

Particle bouncing off the blue wall.

Problem description:

The Tutorial3_BouncingBallElastic.cpp simulates a particle bouncing off a wall. It is assumed the collistion between the particle and the wall is elastic by implying that the restitution coefficient is unity. It means that the particle velocity before and after collision remains the same. Additionally, we will learn how to add a wall over which the ball bounces.

Headers:

The header InfiniteWall.h is included.

#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
#include <Walls/InfiniteWall.h>

Before the main class:

The class InfiniteWall is included in the inherited class Tutorial3:

class Tutorial3 : 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, getZMin()));
        wallHandler.copyAndAddObject(w0);
    }
    
};

Where the object "w0" is an instance of the class InfiniteWall. Then, the function "copyAndAddObject" coppy and add the object to its corresponding handler. The above set of statements, create and place the wall at \( Z_{min} \).
Note: Don’t forget to include the InfiniteWall.h header, as shown in the header section. In some sense, creation and addition of a wall is similar to creation and addition of a particle.

T4: Bouncing ball with dissipation (inelastic)

Particle bouncing off the blue wall with restitution coefficient.

Problem description:

In Tutorial4_BouncingBallInelastic.cpp, the difference between an elastic and inelastic collision between a particle and a wall is illustrated. The only difference between Tutorial3_BouncingBallElastic.cpp and Tutorial4_BouncingBallInelastic.cpp is the value of the restitution coefficient:

double rc = 0.88; // restitution coefficient 

See Bouncing ball - inelastic (code) for more details.

T5: Elastic collision (2 particles)

Particles moving towards each other.

Problem description:
So far, in the above tutorials, we have seen how a particle and a wall interact during a collision. In this tutorial, we illustrate how two particles interact using Tutorial5_ElasticCollision.cpp. For this purpose, we need two particles. The particles may or may not be of the same species type. But, here we shall assume they are of same species and same size.

Before the main function

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
    }
    
};

On comparison between the above class and class Tutorial1, we see how an additional particle is added. Two particles are created, and positioned oppositely apart at a certain distance between them. Both the particles, have a certain velocity directing them towards each other.

T6: Elastic collisions with periodic boundaries

(a) Illustrates the idea behind periodic boundaries, particle exiting boundary b2 re-enters through boundary b1 (b) Illustrates the problem setup.

Problem description:

In order to have multiple collisions, Tutorial5_ElasticCollision.cpp is fitted with periodic boundaries in X.

Headers:

#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
#include <Boundaries/PeriodicBoundary.h>

The header PeriodicBoundary.h is included.

Before the main function:

class Tutorial6 : public Mercury3D
{
public:
    
    void setupInitialConditions() override {
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setRadius(0.005);//particle-1
        p0.setPosition(Vec3D(0.25 * 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);
    }
    
};

In the Class Tutorial6
(i) we create two particles of same type and different sizes. (ii) we setup periodic boundaries in X-direction as

        PeriodicBoundary b0;
        b0.set(Vec3D(1.0, 0.0, 0.0), getXMin(), getXMax());
        boundaryHandler.copyAndAddObject(b0);

Where "b0" is an instance of the class PeriodicBoundary. Then, the object is copied and added to its corresponding handler: boundaryHandler.

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

Particle motion in a box (blue and black denote the walls).

Problem description:

In previous tutorials, we have seen how a particle interacts with a wall and itself. In this tutorial, we will design boxes of different shapes by using more than one wall. As an example, in absence of gravity, we will simulate a particle moving in a two dimensional square shaped box. We consider two dimensions only for the sake of simplicity. The same idea was introduced in Tutorial3_BouncingBallElastic.cpp.

Before the main function:

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);
    }
    
};

In this class, we setup a 2D square shaped box or a polygon by adding more walls as shown above. In total, we have 4 walls forming our box within which the particle will traverse. Note: As we simulate in 2D, no walls are set in z-direction.

Main function:

As our simulation is two dimensional, we set the system dimensions as 2

problem.setSystemDimensions(2);

The complete code for the above problem description can be found in Tutorial7_ParticleIn2DBox.cpp.

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

Problem description:

We extend the problem setup of Tutorial 7, by adding a rectangular block as shown in the above figure. To create this block of wall or obstacle, we will use the Class IntersectionOfWalls. Before we go ahead it is advised to know the difference between an infinite wall and finite wall, see Different types of walls. As an example, we create an obstacle using a set of finite walls and place it within the box created using a set of infinite walls. See the above figure.

Headers:

#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
#include <Walls/InfiniteWall.h>
#include <Walls/IntersectionOfWalls.h>

The header IntersectionOfWalls.h is included.

Before the main function:

To add intersection walls, we use:

        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);

where "w1" is an instance of the class IntersectionOfWalls. This surface (wall) will have the properties of the index 0: "getObject(0)". Then, the object is coppied and added to its corresponding handler: wallHandler.

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

Problem description:

the motion of three particles over an inclined plane is presented in Tutorial9_InclinedPlane.cpp. Here, the particles are fitted with different features.

Headers:

#include <Species/LinearViscoelasticFrictionSpecies.h>
#include <Mercury3D.h>
#include <Walls/InfiniteWall.h>

For Tutorial9_InclinedPlane.cpp, the header which implements the elastic properties and contact force modelc has been changed. In order to include the rolling and sliding effects, LinearViscoelasticFrictionSpecies.h is used.

Before the main function:

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.01 * 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.21 * 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.41 * 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);
    }
};

Three particles are fitted with radius, initial positions and initial velocities.Then, all particles are coppied and added to their corresponding handler.
It is defined the object 'w0' as an instance of the class InfiniteWall. For this tutorial, the surface will have the same properties of "object0". It means the surface and the first particle will have the same elastic properties.

Main function:

int main(int argc, char* argv[])
{
    
    // Problem setup
    Tutorial9 problem;//instantiate an object of class Tutorial 9
    
    double angle = constants::pi / 180.0 * 20.0;
    
    problem.setName("Tutorial9");
    problem.setSystemDimensions(3);
    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.5);
 
    // The normal spring stiffness and normal dissipation is computed and set as
    // For collision time tc=0.005 and restitution coefficeint rc=0.88,
 
    //Properties index 0
    LinearViscoelasticFrictionSpecies species0;
 
    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);
 
    //Properties index 1
    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);
 
    //Combination of properties index 0 and index 1
    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);
 
    //Properties index 2
    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);
 
    //Combination of properties index 0 and index 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);
 
    //Output
    problem.setSaveCount(10);
    problem.dataFile.setFileType(FileType::ONE_FILE);
    problem.restartFile.setFileType(FileType::ONE_FILE);
    problem.fStatFile.setFileType(FileType::ONE_FILE);
    problem.eneFile.setFileType(FileType::NO_FILE);
    
    problem.setXBallsAdditionalArguments("-solidf -v0 -s 8");
    
    problem.setTimeStep(0.005 / 50.0);
    problem.solve(argc, argv);
    
    return 0;
}

The main function of Tutorial9_InclinedPlane.cpp presents how the elastic features of particles are defined. Also, a new function is added to combine features of different indexes: getMixedObject.

T11: Axisymmetric walls inside a cylindrical domain

Problem description:

This tutorial explains how to include axisymmetric walls inside a cylindrical domain. Also, the user will find an useful implementation of how to create polydisperse particle size distribution and perform internal validation tests. The full code can be found in Tutorial11_AxisymmetricWalls.cpp

Headers:

To work with axisymmetric walls, the user needs to add the following headers:

#include "Mercury3D.h"
#include "Walls/IntersectionOfWalls.h"
#include "Walls/AxisymmetricIntersectionOfWalls.h"
#include "Species/LinearViscoelasticSpecies.h"

Before the main function:

The class Tutorial 11 contains the constructor Tutorial 11

    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);
        setParticlesWriteVTK(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);
    }

This constructor is an instance of the main class. Here it's defined the global output and non-changeable parameters.

In this tutorial, it's overrided the member function setupInitialConditions(), which doesn't return any value (defined then as void).

    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));
        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);
        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);
                }
            }
        }
    }

First, it is added the cylinder, then the prism shape or cone that represents the outer part. A flat surface is included in the middle to stop the downforward flow until a specific time. After adding the walls, the species of the particles are included. All particles have the same properties and they are randomly distributed.

An additional member function is included.

    //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()));
        }
    }

This overide function moves the flat surface to the base domain. It happens at a defined time, thus the particles can continue to flow.

Main function:

The main function of this tutorial includes the setup parameters and the tests for a valid collision time and restitution coefficient.

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;
}

T12: Creating objects by using Prisms

Problem description:

This tutorial explains how to define objects in space with the Prism function. Also, the user will find that you can presecribe the location and motion of the object. The full code can be found in Tutorial12_PrismWalls.cpp

Headers:

To work with prismatic walls in this tutorial, the user needs to add the following headers:

#include <Species/LinearViscoelasticSpecies.h>
#include <Mercury3D.h>
#include <Walls/InfiniteWall.h>
#include <Walls/IntersectionOfWalls.h>

Before the main function:

The walls and Prism object are all predefined in the initial conditions.

    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);
        
    }

The first part in the the initial conditions is defining the particle.

        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);

After defining the particle that outer walls are generated by using four finite walls as seen in previous tutorials

To creat the prism we have to do the following: First one creates an instance of IntersectionOfWalls "w1". Next the species are assigned. This line is followed by defining a 3D vector "Points" with space for four vertices. After that the location of the four vertices are defined. Then, by using w1.createPrism(Points) the prism is formed centered at the origin.

        // 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);

Now the object has been formed you can prescribe several properties of the object such as the angular velocity or translational velocity (using setVelocity). Notice that when using setAngularVelocity the object rotates around the origin, so if an object is located at the origin it will rotate around its own axis while located elsewhere it will rotate around the origin. It is essential that you set velocities and positions in the correct order such that you get the desired motion.

        w1.setAngularVelocity(Vec3D(0.0,0.0,1.0));

        w1.setVelocity(Vec3D(0.0,0.0,0.0));

        w1.setPosition(Vec3D(0.5 * getXMax(),0.5 * getYMax(),0.0));

Main:

In the main function the problem is defined in the problem instance, an instance of class Tutorial 12.

    // Problem setup
    Tutorial12 problem; // instantiate an object of class Tutorial 12
    
    problem.setName("Tutorial12");
    problem.setSystemDimensions(2);
    problem.setGravity(Vec3D(0.0, 0.0, 0.0));
    problem.setXMax(0.5);
    problem.setYMax(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 coefficeint 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.setXBallsAdditionalArguments("-solidf -v0 -s .85");
    
    problem.setTimeStep(.005 / 50.0); // (collision time)/50.0

At the end of the main function we call the solve function

    problem.solve(argc, argv);

After solving the problem you can visualize the results and see the effect of the motion of the object on the trajectory of the particle.